.packageName <- "supclust"
## The function implementing the sign-flip for Pelora - empirical covariance
sign.change <- function(x,y)
  {
    if(is.null(dx <- dim(x))) stop("`x' must be a numeric matrix")
    signs <- sign(apply(x, 2, cov, y))
    list(x.new = x * rep(signs, each = dx[1]), signs = signs)
  }


## Computing the coefficients for penalized logistic regression
ridge.coef <- function(x,y,lambda)
  {
    X       <- cbind(rep(1,length(y)),x)
    pnlty   <- diag(apply(X,2,var)*lambda*nrow(X),nr=ncol(X))
    th      <- c(log(mean(y)/(1-mean(y))),rep(0,ncol(X)-1))

    p       <- (1/(1+exp(-(X%*%th))))
    W       <- diag(drop(p*(1-p)))
    th      <- solve((t(X)%*%W%*%X)+pnlty, (t(X)%*%(y-p))+(t(X)%*%W%*%X%*%th))

    p       <- (1/(1+exp(-(X%*%th))))
    W       <- diag(drop(p*(1-p)))
    th      <- solve((t(X)%*%W%*%X)+pnlty, (t(X)%*%(y-p))+(t(X)%*%W%*%X%*%th))
    drop(th)
  }


## The function for the standardization of genes
standardize.genes <- function(exmat)
  {
    means <- apply(exmat,2,mean)
    sdevs <- apply(exmat,2,sd)
    for (i in 1:(dim(exmat)[2]))
      {
        exmat[,i] <- (exmat[,i]-means[i])/sdevs[i]
      }
    list(x=exmat, means=means, sdevs=sdevs)
  }
## Pelora - a supervised algorithm for grouping predictor variables
pelora <- function(x, y, u = NULL, noc = 10, lambda = 1/32, flip = "pm",
                   standardize = TRUE, trace = 1)
  {
    ## Bedeutung der Input-Variablen
    ## -----------------------------
    ## x            Expressionmatrix im Format (n x p)
    ## y            Binrer Responsevektor, codiert mit 0 und 1
    ## u            Klinische Variablen im Format (n x m)
    ## noc          Anzahl Variablen, die ins Modell aufgenommen werden
    ## lambda       Reskalierter Penalty-Parameter, sollte in [0,1] sein
    ## flip         Sign-Flipping in der x-Matrix: stndig, einmalig, nicht
    ## standardize  Variablen-Standardisierung in der x-Matrix: ja oder nein?
    ## trace        Run-Control Output: 0 nichts, 1 moderat, 2 ausfhrlich


    ## Check input
    if(!is.matrix(x) || !is.numeric(x))
        stop("`x' must be a numeric matrix (e.g. gene expressions)")
    n  <- nrow(x)
    if(!is.numeric(y) || length(y) != n || any(y != 0 & y != 1))
        stop(paste("`y' must be a numeric vector of length n = ", nrow(x),
                   "with only 0/1 entries"))
    yvals <- 0:1 # the y-values, aka "class labels"

    ## Sign-Flip
    if (flip=="pm")
      {
        x      <- cbind(x, -x)
        signs  <- c(rep(0,ncol(x)))
      }
    if (flip=="cor")
      {
        sgnChg <- sign.change(x,y)
        x      <- sgnChg$x
        signs  <- sgnChg$signs
      }
    if (flip=="none")
      {
        signs  <- NULL
      }

    ## Standardisierung der x-Variablen
    if (standardize == TRUE)
      {
        stndz  <- standardize.genes(x)
        x      <- stndz$x
        means  <- stndz$means
        sdevs  <- stndz$sdevs
      }
    else
      {
        means <- NULL
        sdevs <- NULL
      }

    ## Setting up the clinical variables
    if((have.u <- !is.null(u))) { ## with or without "clinical variables"?
        if(!is.matrix(u) || !is.numeric(u))
            stop("`u' must be a numeric matrix (e.g. clinical variables)")
        m <- ncol(u)
        E <- cbind(x, u)
    } else {
        m <- 0
        E <- x
    }

    ## Initialisierung
    g         <- ncol(E) # if flip %in% c("cor","none") = px + m
    p         <- noc+1
    px        <- ncol(x)
    X         <- matrix(0, n, p)
    X[,1]     <- 1
    P         <- matrix(0, p, p) ##= diag(apply(X,2,var)) * lambda*n
    theta0    <- log(mean(y)/(1-mean(y)))
    theta     <- c(theta0, rep(0, noc))
    prob      <- rep(1/(1 + exp(-theta0)), n)## == 1/(1 + exp(-X %*% theta))
    W         <- diag(prob*(1-prob))

    ## Earlier options that are now held fixed:
    penflag   <- 0 ## for the standard penalty, !=0 for the nonstandard penalty
    blockflag <- 1 ## for not allowing a gene to enter a cluster more than once
    critflag  <- 0 ## for the log-likelihood criterion, 1 for the L2 criterion
    valiflag  <- 1 ## for validation of genes via backdeletion, 0 to do without

    lSize     <- 2*g*p
    
    ## Aufruf der C-Funktion
    res <- .C("R_clusterer",
              E = 	as.double(E),
              X = 	as.double(X),
              W = 	as.double(W),
              P = 	as.double(P),
              y = 	as.double(y),
              prob = 	as.double(prob),
              theta =   as.double(theta),
              lambda=   as.double(lambda*n),
              n = 	as.integer(n),
              g = 	as.integer(g),
              m = 	as.integer(m),
              p = 	as.integer(p),
              critflag = as.integer(critflag),
              penflag  = as.integer(penflag),
              blockflag= as.integer(blockflag),
              valiflag = as.integer(valiflag),
              traceflag= as.integer(trace),
              ## Output:
              genliste = 	integer(lSize),
              kriterium = 	double (lSize),
              DUP = FALSE,
              PACKAGE = "supclust")[c("genliste", "kriterium")]

    ## Auswertung, bilden der Liste mit Genen und Mittelwerten
    i         <- 1
    genliste  <- res$genliste
    kriterium <- res$kriterium
    genes     <- integer()
    cCrit     <- double()
    alle      <- list()
    crit      <- list()
    while (genliste[i] != 0) {

        if (genliste[i] == -1) {
            i <- i+1
            genes <- c(genes, genliste [i])
            cCrit <- c(cCrit, kriterium[i])
            i <- i+1
        }
        if (genliste[i] == -2) {
            kickout <- genliste[i+1]
            iKick <- which(genes == kickout)
            genes <- genes[-iKick]
            cCrit <- cCrit[-iKick]
            i <- i+2
        }
        if (genliste[i] == -4) { # cluster end
            alle    <- c(alle, list(genes))
            crit    <- c(crit, list(cCrit))
            genes <- integer()
            cCrit <- double()
            i <- i+1
        }
    }
    noclu <- length(alle)

    var.type <- factor(sapply(alle, function(j) any(j > px)),
                       levels = c(FALSE, TRUE),
                       labels = c("Cluster", "Clinical"))
    values <- sapply(alle, function(ids) rowMeans(E[, ids, drop=FALSE]))

    ## FIXME: The samp.names should be attached to `values' / kept
    ## -----  where available;  values should get "predictor names" !!
    dnE <- dimnames(E)
    if(is.null(sNames <- dnE[[1]]))
        sNames <- names(y)
    if(is.null(sNames)) sNames <- as.character(1:n)

    ## Output
    out <- list(genes = alle, values = values, y = y, yvals = yvals,
                lambda = lambda, noc = noc, px = px, flip = flip,
                var.type = var.type, crit = crit, signs = signs,
                samp.names = sNames, gene.names = dnE[[2]], call = match.call())
    class(out) <- "pelora"
    out
  } 


## Returns the fitted values
fitted.pelora <- function(object, ...)
  {
    ## Fitted values
    out <- object$values
    if (is.null(vNames <- colnames(out))) 
      vNames      <- paste("Predictor", 1:object$noc)  
    dimnames(out) <- list(object$samp.names, vNames)
    out
  }


## A short overview of what has been found by Pelora
print.pelora <- function(x, digits = getOption("digits"), details = FALSE, ...)
  {
    ## Preliminaries
    noc      <- x$noc
    fin.crit <- sapply(1:noc, function(i) {cr <- x$crit[[i]]; cr[length(cr)]})
    nGenes   <- unlist(lapply(x$genes,length))
    smDig    <- max(2, digits - 4)
    cCrit    <- format(round(fin.crit, smDig), nsmall = smDig, digits = digits)
    cG       <- format(nGenes)
    cI       <- format(1:noc)

    ## General overview
    cat("\nPelora called with lambda = ", x$lambda, ",", sep="")
    isClust <- x$var.type == "Cluster"
    if (allClust <- all(isClust)) ##  predictors are clusters only
      cat(" ", noc, " cluster", if(noc > 1)"s"," fitted\n\n", sep="")
    else ## predictors are both clusters *and* clinical variables
      cat("\n", (nC <- sum(isClust)), " cluster", if(nC>1)"s", " and ",
          noc - nC, " clinical variable", if ((noc-nC)>1)"s",
          " fitted\n\n", sep="")

    ## Information about each of the predictors
    for (i in 1:noc)
      {
	gic <- x$genes[[i]]
        hic <- x$genes[[i]]
        if (!is.null(x$signs) & any(x$signs)==0) 
          hic[hic>(x$px/2)] <- hic[hic>(x$px/2)]-(x$px/2)
	if (isClust[i])
          {
	    ## the i-th predictor is a gene cluster
	    ng <- length(gic) # == nGenes[i]
	    if(allClust)
              cat("Cluster", cI[i], ": Contains ")
	    else
              cat("Predictor", cI[i], ": Cluster with ")
	    cat(cG[i], " gene", if(ng != 1)"s",
		", final criterion ", cCrit[i],"\n", sep="")

	    if(details)
              {
		## printing all the genes
		cj  <- format(1:ng)
		cgi <- format(hic)
		for (j in 1:ng)
                  {
                    cat("Entry", cj[j], ": Gene", cgi[j])
                    if (!is.null(x$signs))
                      {
                        if (x$signs[gic[j]] == 0)
                          {
                            cat(if (gic[j]>(x$px/2))
                                " (flipped)" else "          ")
                          }
                        else
                          {           
                            cat(if (x$signs[gic[j]] == -1)
                                " (flipped)" else "          ")
                          }
                        if (!is.null(x$gene.names))
                          cat(" : Name", x$gene.names[gic[j]])
                        cat("\n")
                      }
                  }
              }
          }
	else
          { ## x$var.type[i] == "Clinical"
	    cat("Predictor ", cI[i], " : Clinical variable ",
                gic[1]-x$px,
		if (!is.null(x$gene.names))
		paste(" named `", x$gene.names[gic[1]], "'"),
		", final criterion ", cCrit[i], "\n", sep="")
          }
	if(details) cat("\n")
      } ## for(j )
    cat("\n")
    invisible(x)
  }

## Yields printed output about the clustering in more detail
summary.pelora <- function(object, digits = getOption("digits"), ...)
  {
    print(object, details = TRUE, digits = digits, ...)
  }


## Plots a 2-dimensional projection of the first 2 Pelora-clusters
plot.pelora <- function(x, main = "2-Dimensional Projection Pelora's output",
                        xlab = NULL, ylab = NULL, col = seq(x$yvals), ...)
  {
    if(x$noc <= 1)
      {
        if(is.null(xlab))
          xlab <- if (x$var.type[1] == "Cluster")
            "Mean Expression of Pelora's Cluster 1"
          else "Expression of a Clinical Variable"

        if(is.null(ylab))
          ylab <- "Class label"
      
        plot(x$values[,1], x$y, type="n", xlab=xlab, ylab=ylab, main=main, ...)
        text(x$values[,1], x$y, x$y, col = col[1 + x$y])
        return()
      }
    
    if(is.null(xlab))
    xlab <-
      if (x$var.type[1] == "Cluster")
        "Mean Expression of Pelora's Cluster 1"
      else "Expression of a Clinical Variable"

    if(is.null(ylab))
    ylab <-
        if (x$var.type[2] == "Cluster")
            "Mean Expression of Pelora's Cluster 2"
        else "Expression of a Clinical Variable"

    plot(x$values[,1], x$values[,2], type = "n",
         xlab = xlab, ylab = ylab, main = main, ...)
    text(x$values[,1], x$values[,2], x$y, col = col[1 + x$y])
    invisible()
}


## Returns the coefficients of the penalized logistic regression classifier
coef.pelora <- function(object, ...)
  {
    clusnames  <- character(object$noc)
    for (i in 1:object$noc) clusnames[i] <- paste("Predictor", i)
    out        <- ridge.coef(fitted(object), object$y, object$lambda)
    names(out) <- c("Intercept", clusnames)
    out
  }


## Predictions with Pelora's clusters
predict.pelora <- function(object, newdata = NULL, newclin = NULL,
                           type = c("fitted", "probs", "class"),
                           noc = object$noc, ...)
  {
    ## Checking the input
    type <- match.arg(type)

    ## Return fitted values, probabilities or class labels for the training data
    if (is.null(newdata))
      {
        X   <- cbind(rep(1,length(object$y)),fitted(object))
        
        if (length(noc)==1 && noc>object$noc)
          stop("You cannot predict with more predictors than you have fitted")
        
        if (length(noc)==1 && noc==object$noc)
          {
            prvec           <- matrix(1/(1+exp(-c(X%*%coef(object)))), ncol = 1)
            dimnames(prvec) <- list(object$samp.names, paste(noc, "Predictors"))
            return(switch(type,
                          fitted = fitted(object),
                          probs  = prvec,
                          class  = (prvec > 0.5)*1))
          }

        if (length(noc)==1 && noc<object$noc)
          {
            koeff <- ridge.coef(object$val[, 1:noc], object$y, object$lambda)
            prvec <- matrix(c((1/(1+exp(-(X[, 1:(noc+1)]%*%koeff))))), ncol = 1)
            dimnames(prvec) <- list(object$samp.names, paste(noc, "Predictors"))
            return(switch(type,
                          fitted = fitted(object),
                          probs  = prvec,
                          class  = (prvec>0.5)*1))
          }

        if (length(noc)>1 & max(noc)>object$noc)
          stop("You cannot predict with more predictors than you have fitted")
        
        if (length(noc)>1 & max(noc)<=object$noc)
          {
            prmt <- NULL
            for (i in 1:length(noc)) {
              koeff <- ridge.coef(object$val[,1:(noc[i])],object$y,object$lamb)
              prvec <- matrix(c((1/(1+exp(-(X[,1:(noc[i]+1)]%*%koeff))))),nc=1)
              dimnames(prvec) <- list(object$samp.n,paste(noc[i],"Predictors"))
              prmt  <- cbind(prmt, prvec)
            }
            return(switch(type,
                          fitted = fitted(object)[, noc, drop = FALSE],
                          probs  = prmt,
                          class  = (prmt>0.5)*1))
          }
      }

    ## Returning fitted values, probabilities or class labels for test data
    else
      {
        ## Customizing newdata according to the choice of the flipping method
        if (object$flip=="pm") newdata <- cbind(newdata, -newdata)

        ## Check if new clinical variables are provided too
        if (is.null(newclin) && any(object$var.type=="Clinical"))
          stop("You also need to provide new clinical variables")
                
        ## Check the dimensions of the new data
        if (ncol(newdata)!=object$px)
          stop(paste("The new data need to have the same number of",
                     "predictor variables as the training data"))

        ## Flip the signs of the new data
        if (object$flip=="cor") newdata <- t(t(newdata)*object$signs)

        ## Standardization of the new data
        if (!is.null(object$means))
          {
            for (i in 1:(ncol(newdata)))
              {
                newdata[,i] <- (newdata[,i]-object$means[i])/object$sdevs[i]
              }
          }

        ## Merge expression data and clinical variables (if available)
        if (!is.null(newclin)) newdata <- cbind(newdata, newclin)
         
        ## Determine the fitted values
        Xt <- cbind(rep(1,nrow(newdata)))
        for (j in 1:object$noc)
          {
            Xt <- cbind(Xt, rowMeans(newdata[,object$genes[[j]], drop = FALSE]))
          }
        
        ## Naming the predictors
        sampnames <- 1:nrow(newdata)
        clusnames <- character(object$noc)
        for(i in 1:object$noc)      clusnames[i] <- paste("Predictor", i)
        dimnames(Xt) <- list(sampnames, c("Intercept", clusnames))

        if (length(noc)==1 && noc>object$noc)
          stop("You cannot predict with more predictors than you have fitted")
        
        if (length(noc)==1 && noc==object$noc)
          {
            prvec <- matrix(1/(1+exp(-c(Xt%*%coef(object)))), ncol=1)
            dimnames(prvec) <- list(1:nrow(newdata), paste(noc, "Predictors"))
            return(switch(type,
                          fitted = Xt[,2:ncol(Xt)],
                          probs  = prvec,
                          class  = (prvec > 0.5)*1))
          }

        if (length(noc)==1 && noc<object$noc)
          {
            koeff <- ridge.coef(object$val[,1:noc],object$y,object$lambda)
            prvec <- matrix(c((1/(1+exp(-(Xt[,1:(noc+1)]%*%koeff))))), ncol=1)
            dimnames(prvec) <- list(1:nrow(newdata), paste(noc, "Predictors"))
            return(switch(type,
                          fitted = Xt[,2:ncol(Xt)],
                          probs  = prvec,
                          class  = (prvec>0.5)*1))
          }

        if (length(noc)>1 && max(noc)>object$noc)
          stop("You cannot predict with more predictors than you have fitted")
        
        if (length(noc)>1 && max(noc)<=object$noc)
          {
            prmt <- NULL
            for (i in 1:length(noc)) {
              koeff <- ridge.coef(object$val[,1:(noc[i])],object$y,object$lamb)
              prvec <- matrix(c((1/(1+exp(-(Xt[,1:(noc[i]+1)]%*%koeff))))),nc=1)
              dimnames(prvec)<- list(1:nrow(newdata),paste(noc[i],"Predictors"))
              prmt  <- cbind(prmt, prvec)
            }
            return(switch(type,
                          fitted = Xt[,2:ncol(Xt)],
                          probs  = prmt,
                          class  = (prmt>0.5)*1))
          }
      }
  } 













## The function implementing the cleaning stage for Wilma
back.search <- function(genes, x, y, verbose = FALSE)
{
    reduced <- TRUE
    while (reduced && (N <- length(genes)) >= 2)
    {
	reduced <- FALSE
	if (score(mx <- rowMeans(x[, genes]), y) == 0)
	{
            margins <- sapply(seq(along=genes), function(i)
                              margin(rowMeans(x[, genes[-i], drop=FALSE]), y))
	    if (max(margins) > margin(mx, y)) {
                imax <- which.max(margins)
		if(verbose)
		    cat("Gen: ", genes[imax],
			"Score:	  0",
			"Margin: ", round(max(margins),3),"\n")
		genes	   <- genes[-imax]
		reduced	   <- TRUE
	    }
	}
	
	if ((sc <- score(mx <- rowMeans(x[, genes]), y)) > 0)
	{
	    old.score <- sc
            scores	<- numeric(N)
	    for (i in 1:length(genes))
		scores[i]  <- score(rowMeans(x[, genes[-i], drop=FALSE]), y)

            minsc <- min(scores)
	    indices <- which(scores == minsc)# all of them
            mgs <- sapply(indices, function(ii)
                          margin(rowMeans(x[, genes[- ii], drop=FALSE]), y))

	    if (minsc < old.score || max(mgs) > margin(mx, y)) {
                imax <- which.max(mgs)
		if(verbose)
		    cat("Gen: ", indices[imax],
			"Score: ", minsc,
			"Margin: ", round(mgs[imax],3),"\n")
		genes	   <- genes[-indices[imax]]
		reduced	   <- TRUE
              }
          } 
      }
    genes
  }


## The margin function for Wilma
margin <- function(x, resp)
  {
    .C("R_margin",
       as.double(x[order(resp)]),
       as.integer(sum(resp==0)),
       as.integer(sum(resp==1)),
       re = double(1),
       PACKAGE = "supclust")$re
  }


## The score function for Wilma
score <- function(x, resp)
  {
    .C("R_score",
       as.double(x),
       as.integer(resp),
       as.integer(length(x)),
       re = double(1),
       PACKAGE = "supclust")$re
  }


## The function implementing the sign-flip for Wilma -- using score()
sign.flip  <- function(x, y)
  {
    y <- as.integer(y)
    O <- as.integer(0)
    if(any(y < O | y > 1:1)) stop("`y' must be 0 or 1!")
    n <- length(y)
    scores <- apply(x, 2, score, y)
    n0 <- sum(y == O)
    middle <- n0 * (n - n0) / 2
    lrg  <- scores > middle
    x[,lrg] <- -x[,lrg]
    list(flipped.matrix = x, signs = ((-2)*lrg)+1)
  }


## Classification with Wilma's clusters - nearest neighbor rule
nnr <- function(xlearn, xtest, ylearn)
  {
    as.numeric(knn(xlearn, xtest, ylearn))-1
  }


## Classification with Wilma's clusters - diagonal linear discriminant analysis
dlda <- function (xlearn, xtest, ylearn)
  {
    ## Definition of variables    
    n      <- nrow(xlearn)
    p      <- ncol(xlearn)
    nk     <- rep(0, max(ylearn) - min(ylearn) + 1)
    K      <- length(nk)
    m      <- matrix(0, K, p)
    v      <- matrix(0, K, p)
    disc   <- matrix(0, nrow(xtest), K)

    ## Computing mean and variances
    for (k in (1:K))
      {
        which   <- (ylearn == k + min(ylearn) - 1)
        nk[k]   <- sum(which)
        m[k, ]  <- apply(xlearn[which, , drop = FALSE], 2, mean)
        v[k,]   <- apply(xlearn[which, , drop = FALSE], 2, var)
      }

    ## Computing the pooled variance
    vp <- apply(v, 2, function(z) sum((nk - 1) * z)/(n - sum(nk!=0)))

    ## Computing the discriminant function
    for (k in (1:K))
      {
        disc[,k] <- apply(xtest, 1, function(z) sum((z-m[k,])^2*(1/vp)))
      }

    ## Prediction and output
    pred <- apply(disc, 1, function(z) (min(ylearn):max(ylearn))[order(z)[1]])
    pred
   }


## Classification with Wilma's clusters - logistic regression
logreg <- function(xlearn, xtest, ylearn)
  {
    xvalues    <- xlearn
    op         <- options(warn = -1)
    model      <- glm(ylearn~., data = data.frame(xvalues), family = binomial)
    xvalues    <- xtest
    predic     <- predict(model,new = data.frame(xvalues), type = "response")
    options(op)
    as.numeric(predic >= 0.5)
  }


## Classification with Wilma's clusters - aggregated trees
aggtrees <- function(xlearn, xtest, ylearn)
  {
    noc     <- ncol(xtest) 
    predic  <- matrix(0, nrow(xtest), noc)
    for (j in 1:noc)
      {
        xvalues    <- xlearn
        model      <- rpart(ylearn ~ ., data = data.frame(xvalues))
        xvalues    <- xtest
        predic[,j] <- predict(model, newdata = data.frame(xvalues))
      } 
    final <- rowSums(predic)
    for (j in 1:nrow(xtest))
      {
        if (final[j] == (noc/2)) final[j] <- final[j]+(predic[j,1]-0.5)
      }
    as.numeric(final > (round((noc-0.1)/2)))
  }

## Wilma - an algorithm for supervised grouping of predictor variables
wilma <- function(x, y, noc, genes = NULL, flip = TRUE,
                  once.per.clust = FALSE, trace = 0)
{
    ## Checking the input
    y <- as.integer(y)
    if(any(y < 0 | y > 1))
        stop("Labels y have to be 0 or 1!")
    n <- length(y)
    if(!is.matrix(x) || !is.numeric(x) || nrow(x) != n)
        stop("`x' must be a numeric matrix with `n' (= length(y)) observations")
    if((noc <- as.integer(noc)) < 1)
        stop("`noc' must be a positive integer")
    if((trace <- as.integer(trace)) < 0)
        stop("`trace' must be integer >= 0 (or logical)")
    ## C output (Cverb > 0) only for trace >= 2 :
    Cverb <- as.integer(if(trace) Cverb <- trace - 1 else 0)

    ## Customizing the problem and sign-flipping
    iy    <- sort.list(y)# i.e., first the 0's, then the 1's
    io    <- order(iy)   
    y     <- y[iy]
    x     <- x[iy,]
    signs <- NULL
    if (flip)
      {
        res   <- sign.flip(x,y)
        x     <- res$flipped.matrix
        signs <- res$signs
      }

    ## Definitions
    n1 <- sum(y == 0)
    n2 <- n - n1 ## = sum(y == 1)
    p  <- ncol(x) ## = length(x)/n


    ## Looking for (optional) starting genes
    used <- rep(FALSE, p)
    if(!is.null(genes)) {
        if(!is.list(genes) || length(genes) != noc)
            stop("starting `genes' must be a list of length `noc' (=", noc,")")
        sgl <- unlist(genes)
        ##was sgl <- NULL ; for (i in 1:noc) sgl   <- c(sgl, genes[[i]])
        used[sgl] <- TRUE
    }

    mn.x <- gList <- vector("list", noc)
    steps <- integer(noc)

    ## The loop for supervised clustering
    for (i in 1:noc)
    {
        ## Initial value if (optional) starting clusters are provided
        size <- length(gic <- genes[[i]])
        clMean <- rep(0,n)

        if(trace) cat("\n", "\nCluster ", i, "\n----------\n", sep="")

        ## Forward search and cleaning stages
        nFwd <- 1
        repeat { ## Search forward and backward once :

            if (size > 0)
                clMean <- rowMeans(x[, gic, drop=FALSE])

            if(trace && nFwd > 1) {
                cat("used[]:", which(as.logical(used)),"\n")
                if(trace >= 2 && size > 0)
                    cat("entry gic[]:", gic, "\n")
            }

            res <- .C("R_multicluster",
                      as.double(x), as.integer(y),# 2
                      as.integer(n), as.integer(n1), as.integer(n2),# 5
                      as.integer(p),# 6
                      used = as.integer(used),# 7
                      as.double(clMean),
                      glsize = size,# 9
                      gic    = c(as.integer(gic), integer(p)),
                      scores = integer(size+p),
                      margins=  double(size+p),
                      as.logical(once.per.clust),
                      Cverb,
                      PACKAGE="supclust")[
                      c("used", "glsize", "gic", "scores", "margins")]

            if(trace) {
                if(nFwd > 1 && trace >= 2 && size > 0)
                    cat("exit gic[1:res$glsize]:", res$gic[1:res$glsize], "\n")
                if(res$glsize > size) {
                    cat("\nAccepted",if(size > 0) "(additionally)","\n")
                    for(j in (size+1):res$glsize)
                        .p.1gen(j,  res$ gic [j],
                                s = res$ scores[j],
                                m = res$margins[j], digits = 3)
                    cat(" gic size changed from ", size," to ",res$glsize,"\n")
                }
                else ## res$glsize == size
                    cat(" gic[] *unchanged*\n")
            }
            gic.old <- gic
            size <- res$glsize
            gic  <- res$gic[1:size]
            used <- res$used

            if(trace) cat("\nEliminating ")
            gic.red  <- back.search(gic, x, y, verbose = trace > 0)
            if(length(gic.red) == size) {
                if(trace)
                    cat(" -- no reduction --> end{repeat} after ",
                        nFwd, if(nFwd == 1)"step" else "steps","\n")
                break
            }

            ## else : *have* reduced

            w <- gic[!(gic %in% gic.red)]
            if(trace) cat(" w= (", paste(w,collapse=','),")", sep="")
            used[w] <- FALSE
            gic <- gic.red
            size <- length(gic)

            if(trace)
                cat(" -- end {repeat} nr. ", nFwd,
                    ". gic enlarged and reduced from ", length(gic.old),
                    " to ", size, "\n", sep="")

            if (nFwd > 1 && length(gic.old) == size && all(gic.old == gic))
                break
            nFwd <- nFwd + 1

        }## end{repeat}

        gList[[i]] <- gic
        steps [i] <- nFwd
        mn.x [[i]] <-
            if(size) sapply(1:size,
                            function(j) rowMeans(x[, gic[1:j], drop=FALSE]))
            else integer(0)

	if(trace)
	    p.1clust(i, gic = gic, x.mean = mn.x[[i]], y = y, nFwd = nFwd)

    }## end for i = 1:noc

    ## Restoring the original order of output and response
    y <- y[io]
    for (i in 1:length(mn.x)) mn.x[[i]] <- mn.x[[i]][io,]

    ## Output
    r <- list(clist = gList,
              steps = steps, y = y, x.means = mn.x, noc = noc, signs=signs)
    class(r) <- "wilma" # cluster List
    r
}

## Short overview of what has been done by Wilma
print.wilma <- function(x, ...)
  {
    ## The number of clusters which were fitted
    noc <- x$noc

    ## Final score and margin values
    fvals <- fitted(x)
    s     <- apply(fvals, 2, score, x$y)
    m     <- apply(fvals, 2, margin, x$y)

    ## Formatting the numbers
    nGenes <- unlist(lapply(x$clist,length))
    cScore <- format(s)
    cMargi <- format(round(m,2))
    cG     <- format(nGenes)
    cI     <- format(1:noc)

    ## General overview
    cat("\nWilma called to fit ", x$noc, " cluster", if(noc>1) "s",
        "\n\n", sep="")

    ## Information for each of the predictors
    for (i in 1:noc)
      {
	gic <- x$clist[[i]]
        ng <- length(gic) 
        cat("Cluster", cI[i], ": Contains ")
        cat(cG[i], " gene", if(ng != 1)"s", ", final score ", cScore[i],
            ", final margin ", cMargi[i], "\n", sep="")
      }
  }


## The next 3 functions provide a detailed overview about Wilma's clustering
summary.wilma <- function(object, ...)
  {
    cat("`Wilma' object: ")
    printClist(object, ...)
  }

## Auxiliary function for summary.wilma()
printClist <- function(x, ...)
  {
    ## Purpose: print (Method) "gene list" objects, called from print.multiclust
    ## ----------------------------------------------------------------------
    ## Author: Martin Maechler, Date: 21 Jun 2003, 21:38

    if(is.null(gic <- x$clist))
        stop("invalid `x' argument")
    noc <- length(gic)
    cat("number of clusters `noc' =", noc,"\n")
    for(i in 1:noc)
        p.1clust(i, gic = gic[[i]],
                 x.mean = x$x.means[[i]], y = x$y, nFwd = x$steps[i])
    invisible(x)
  }

## Another auxiliary function for summary.wilma()
p.1clust <- function(i, gic, x.mean, y, nFwd)
  {
    cat("\nFinal Cluster", i,
        "\n----------------\n\n")
    size <- length(gic)
    if(size >= 1)
        for (j in 1:size) {
            s <- score (x.mean[,j], y)
            m <- margin(x.mean[,j], y)
            .p.1gen(j=NULL, id = gic[j], s, m, digits = 3)
        }
    else cat(" __empty__ \n")
    invisible()
  }

.p.1gen <- function(j, id, s, m, digits = 3)
  {
    cat("Gen",
        if(is.null(j)) "" else paste("[", j, "]",if(j <= 9)" ",sep=""),
        ":", formatC(id, wid=6),
        "  Score:", formatC(s, wid=4),
        "  Margin:",
        formatC(format(c(round(m, digits=digits),0.735))[1], wid=7, dig=3),
        "\n", sep="")
  }



## Returns the fitted values
fitted.wilma <- function(object, ...)
  {
    ## Fitted values
    out <- NULL
    for (i in 1:object$noc)
      {
        out <- cbind(out, (object$x.means[[i]])[,ncol(object$x.means[[i]])])
      }
    
    if (is.null(vNames <- colnames(out))) ## Naming the predictors
        vNames <- paste("Predictor", 1:object$noc)
    dimnames(out) <- list(1:length(object$y), vNames)
    out
}

## 2-dimensional projection of Wilma's output
plot.wilma <- function(x, xlab = NULL, ylab = NULL, col = seq(x$yvals),
                       main = "2-Dimensional Projection of Wilma's Output", ...)
  {
    ## Fitted values
    fvals <- fitted(x)

    ## If only 1 cluster was fitted
    if(x$noc <= 1)
      {
        if(is.null(xlab)) xlab <- "Mean Expression of Wilma's Cluster 1"
        if(is.null(ylab)) ylab <- "Class label"
        plot(fvals[,1], x$y, type = "n", xlab=xlab, ylab=ylab, main=main, ...)
        text(fvals[,1], x$y, x$y, col = col[1 + x$y])
        return()
      }
    
    ## For 2 clusters and more
    if(is.null(xlab)) xlab <- "Mean Expression of Wilma's Cluster 1"
    if(is.null(ylab)) ylab <- "Mean Expression of Wilma's Cluster 2"
    plot(fvals[,1], fvals[,2], type="n", xlab=xlab, ylab=ylab, main=main, ...)
    text(fvals[,1], fvals[,2], x$y, col=col[1 + x$y])
    invisible()
}


predict.wilma <- function(object, newdata = NULL, type = c("fitted", "class"),
                          classifier = c("nnr", "dlda", "logreg", "aggtrees"),
                          noc = object$noc, ...)
  {
    ## Checking the input
    type       <- match.arg(type)
    classifier <- match.arg(classifier)
    
    ## Return fitted values for the training data
    if (is.null(newdata))
      {
        if (length(noc)==1 && noc>object$noc)
            stop("You cannot predict with more predictors than you have fitted")
                 
        if (length(noc)==1 && noc<=object$noc)
            return(fitted(object))
        
        if (length(noc)>1 & max(noc)>object$noc)
            stop("You cannot predict with more predictors than you have fitted")
        
        if (length(noc)>1 & max(noc)<=object$noc)
            return(fitted(object)[, noc, drop = FALSE])
      }

    ## Returning fitted values or class labels for test data
    else
      {
        ## Check the dimensions of the new data
        if (!is.null(object$signs) && ncol(newdata)!=length(object$signs))
          {
            stop(paste("The new data need to have the same number of",
                       "predictor variables as the training data"))
          }

        ## Flip the signs of the new data
        if (!is.null(object$signs))
          {
            newdata <- t(t(newdata)*object$signs)
          }
        
        ## Fitted values for the new data
        xtest <- NULL
        for (j in 1:object$noc)
          {
            xtest <- cbind(xtest, rowMeans(newdata[,object$clist[[j]],
                                                   drop = FALSE]))
          }

        ## Naming the predictors and samples
        sampnames <- 1:nrow(newdata)
        clusnames <- character(object$noc)
        for(i in 1:object$noc)      clusnames[i] <- paste("Predictor", i)
        dimnames(xtest) <- list(sampnames, clusnames)

        ## Return the fitted values for the new data if requested
        if (type == "fitted") return(xtest[,noc, drop = FALSE])

        ## Do 0/1-classification
        if (length(noc)==1 && noc>object$noc)
            stop("You cannot predict with more predictors than you have fitted")
        
        if (length(noc)==1 && noc==object$noc)
          {
            xlearn <- fitted(object)
            return(switch(classifier,
                          nnr      = nnr(xlearn, xtest, object$y),
                          dlda     = dlda(xlearn, xtest, object$y),
                          logreg   = logreg(xlearn, xtest, object$y),
                          aggtrees = aggtrees(xlearn, xtest, object$y)))
          }

        if (length(noc)==1 && noc<object$noc)
          {
            xlearn <- (fitted(object))[, 1:noc, drop = FALSE]
            xtest  <- xtest[, 1:noc, drop = FALSE]
            return(switch(classifier,
                          nnr      = nnr(xlearn, xtest, object$y),
                          dlda     = dlda(xlearn, xtest, object$y),
                          logreg   = logreg(xlearn, xtest, object$y),
                          aggtrees = aggtrees(xlearn, xtest, object$y)))
          }

        if (length(noc)>1 && max(noc)>object$noc)
            stop("You cannot predict with more predictors than you have fitted")
          
        if (length(noc)>1 && max(noc)<=object$noc)
          {
            clmt <- NULL
            for (i in 1:length(noc))
              {
                xl <- (fitted(object))[, 1:(noc[i]), drop = FALSE]
                xt <- xtest[, 1:(noc[i]), drop = FALSE]
                cl <- switch(classifier,
                             nnr      = nnr(xl, xt, object$y),
                             dlda     = dlda(xl, xt, object$y),
                             logreg   = logreg(xl, xt, object$y),
                             aggtrees = aggtrees(xl, xt, object$y))
                cl           <- matrix(cl, ncol=1)
                dimnames(cl) <- list(1:nrow(newdata),paste(noc[i],"Predictors"))
                clmt         <- cbind(clmt, cl)
              }
            return(clmt)
          }
      }
  } 


.First.lib <- function(lib, pkg) {
    stopifnot(require(class)) # for knn()
    stopifnot(require(rpart)) # for rpart()
    library.dynam("supclust", pkg, lib)
}


## for R versions < 1.5
if(paste(R.version$major, R.version$minor, sep=".") < 1.5) {
    ## cheap substitutes :
    rowMeans <- function(x) apply(x, 1, mean) ## used in many places
    colMeans <- function(x) apply(x, 2, mean)
    rowSums  <- function(x) apply(x, 1, sum)
    colSums  <- function(x) apply(x, 2, sum)
}
