
cclust <- function (x, centers, iter.max=100, verbose=FALSE, dist="euclidean",
                    method= "kmeans", rate.method="polynomial", rate.par=NULL)
{
    
  xrows<-dim(x)[1]
  xcols<-dim(x)[2]
  
  perm <- sample(xrows)
  x <- x[perm,]
  
  if (is.matrix(centers))   # initial values are given
      ncenters <- dim(centers)[1]
  else
    {                   # take centers random vectors as initial values
      ncenters <- centers
      centers <- x[rank(runif(xrows))[1:ncenters],]
    }

  dist <- pmatch(dist, c("euclidean", "manhattan"))
            if (is.na(dist)) 
                    stop("invalid distance")
            if (dist == -1) 
                    stop("ambiguous distance")
         
  
  method <- pmatch(method, c("kmeans", "hardcl", "neuralgas"))
            if (is.na(method)) 
                    stop("invalid clustering method")
            if (method == -1) 
                    stop("ambiguous clustering method")
           

 rate.method<- pmatch(rate.method, c("polynomial", "exponentially.decaying"))
            if (is.na(rate.method)) 
                    stop("invalid learning rate method")
            if (rate.method == -1) 
                    stop("ambiguous learning rate method")
           

  if (method==2) {
   if (rate.method==1 && missing (rate.par)) {
     rate.par<-c(1.0,0.0)
   }
   else if (rate.method==2 && missing (rate.par)) {
     rate.par<-c(1e-1,1e-4)
   }
 }
   if (method==3 && missing (rate.par)) {
     rate.par<-c(0.5,0.005,10,0.01)
   }
  
  initcenters <- centers
  dist <- matrix(0,xrows,ncenters)
  pos <- as.factor(1:ncenters)   # necessary for empty clusters
  rownames(centers) <- pos

  iter <- integer(1)
  changes <- integer(iter.max)
  cluster <- integer(xrows)
  clustersize <- integer(ncenters)


  if (method==1){
  retval <- .C("kmeans",
               xrows = as.integer(xrows),
               xcols = as.integer(xcols),
               x = as.double(x),
               ncenters = as.integer(ncenters),
               centers = as.double(centers),
               cluster = as.integer(cluster),
               iter.max = as.integer(iter.max),
               iter = as.integer(iter),
               changes = as.integer(changes),
               clustersize = as.integer(clustersize),
               verbose = as.integer(verbose),
               dist = as.integer(dist)) 
   
 } 
  else if (method==2){
  retval <- .C("hardcl",
               xrows = as.integer(xrows),
               xcols = as.integer(xcols),
               x = as.double(x),
               ncenters = as.integer(ncenters),
               centers = as.double(centers),
               cluster = as.integer(cluster),
               iter.max = as.integer(iter.max),
               iter = as.integer(iter),
               clustersize = as.integer(clustersize),
               verbose = as.integer(verbose),
               dist = as.integer(dist),
               methrate = as.integer(rate.method),
               par=as.double(rate.par))
}
  else if (method==3){
  retval <- .C("neuralgas",
               xrows = as.integer(xrows),
               xcols = as.integer(xcols),
               x = as.double(x),
               ncenters = as.integer(ncenters),
               centers = as.double(centers),
               cluster = as.integer(cluster),
               iter.max = as.integer(iter.max),
               iter = as.integer(iter),
               clustersize = as.integer(clustersize),
               verbose = as.integer(verbose),
               dist = as.integer(dist),
               par=as.double(rate.par))
  
}
  centers <- matrix(retval$centers,
                    ncol=xcols, dimnames=dimnames(initcenters))
  cluster <- retval$cluster + 1
  cluster <- cluster[order(perm)]
  
  if (method==1) {
    methrate <- NA
    par <- NA
  }
  if (method==3){
    methrate <- NA
  }
  
  retval <- list (centers = centers,
                  initcenters = initcenters,
                  ncenters = ncenters,
                  cluster = cluster,
                  size = retval$clustersize,
                  iter = retval$iter - 1,
                  changes = retval$changes,
                  dist = dist,
                  method = method,
                  rate.method=rate.method,
                  rate.par=rate.par,
                  call=match.call())

  class(retval)<-c("cclust")
  return(retval)
}


print.cclust <- function (clobj)
  {
  
    if (!is.null(clobj$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$ncenters, "\n")
    cat("Sizes  of Clusters: ", clobj$size, "\n\n")
    cat("Learning Parameters:",clobj$rate.par,"\n\n")
    
  if (clobj$method==1)
    {
  if (!is.null(clobj$iter))
      {
        if (clobj$iter < length(clobj$changes))
          cat("Algorithm converged after", clobj$iter, "iterations.\n")
        else
          cat("Algorithm did not converge after", clobj$iter, "iterations.\n")
        cat("Changes:", clobj$changes[1:clobj$iter], "\n\n")
      }
    }
 
  }

plot.cclust <- function(clobj, x, centers=TRUE, initcenters=TRUE,
                         color=rainbow(clobj$ncenters),...){
  
  x <- as.matrix(x)
  
  cl <- predict(clobj, 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="X",col=color)
    if(initcenters)
      points(cl$initcenters, pch="+",col=color)
  }
}


predict.cclust <- function(clobj, x){

  xrows<-dim(x)[1]
  xcols<-dim(x)[2]
  ncenters <- clobj$ncenters
  cluster <- integer(xrows)
  clustersize <- integer(ncenters)
  

  if(dim(clobj$centers)[2] != xcols){
    stop("Number of variables in cluster object and x are not the same!")
  }

  
  retval <- .C("assign",
               xrows = as.integer(xrows),
               xcols = as.integer(xcols),
               x = as.double(x),
               ncenters = as.integer(ncenters),
               centers = as.double(clobj$centers),
               cluster = as.integer(cluster),
               clustersize = as.integer(clustersize),
               dist = as.integer(clobj$dist))

  
     

  clobj$initcenters <- NULL
  clobj$iter <- NULL
  clobj$changes <- NULL
  clobj$cluster <- retval$cluster+1
  clobj$size <- retval$clustersize

  return(clobj)
}











.First.lib <- function(lib, pkg) library.dynam("cclust", pkg, lib)
