iterative.train.SOM <- function(data, codebook, radius, 
                                 train.steps, alpha, Coords, H.distance) {
      #========================================
      # data         - the data on which the SOM is to be run
      # codebook     - the codebook (prototype vectors)
      # radius       - the radius (in output space) 
      # train.steps
      # alpha
      # Coords
      # H.distance
      #========================================
   update.step <- 100
      # number of cases in data
   data.len <- dim(data)[1]           #300
   trainlen <- train.steps
   steps <- 1:trainlen
   samples <- sample(1:data.len, update.step)
   radius <- radius^2
   radius[radius == 0] <- 2.2204e-016       # ~ Machine eps
   for (t in 1:trainlen) {
      x <- data[samples[t], ]              # pick one sample vector
      Dx <-  t(t(codebook) - x)             # get the difference from codebook vectors

      bmu <- it.min.dist(x, codebook)$where  # Find the closest
           # For the BMU get an adjustment to the values for the codebook vectors
           # based on the proximity to the BMU in the output space 
      h <- alpha[t] * exp(-H.distance[,bmu]/(2*radius[t])) 
           # update codebook
      codebook <- codebook - h*Dx
   }
   return(codebook)
}

it.min.dist <- function(pt, array) {
   D <- sqrt(apply((t(t(array)-as.vector(pt,mode="numeric")))^2, 1, sum))
   return(list(min=min(D),where=which(min(D)==D)))
}

N.iterative.train.SOM <- function(data, proto.vectors, radius, 
                                 trainlen, alpha, Coords, H.distance, pass) {
      #========================================
      # data         - the data on which the SOM is to be run
      # proto.vectors     - the proto.vectors (prototype vectors)
      # radius       - the radius (in output space) 
      # trainlen
      # alpha
      # Coords
      # H.distance
      #========================================
      # number of cases in data
   data.len <- dim(data)[1]           
       # Normally random
   samples <- seq(1,21,by=5)+pass
   rad <- radius[1]
   radius <- radius^2
   radius[radius == 0] <- 2.2204e-016       # ~ Machine eps
   for (t in 1:trainlen) {
      x <- data[samples[t], ]        # pick one sample vector
      points(x[1],x[2], col="blue", cex=2)
      Dx <-  t(t(proto.vectors) - x)      # get the difference from proto.vectors vectors
      bmu <- it.min.dist(x, proto.vectors)$where  # Find the closest
#      points(proto.vectors[bmu,1], proto.vectors[bmu,2],pch=16, col="green")
                  # For the BMU get an adjustment to the
                  # values for the proto.vectors vectors
                  # based on the proximity to the BMU in
                  # the output space
                  # The H.distance[,bmu] gives the distance on the hex map
                  # from the BMU to all other pv's
      h <- alpha[t] * exp(-H.distance[,bmu]/(2*radius[t]))
      h.Dx <- h*Dx
           # update proto.vectors based on distance
      ind <- order(H.distance[,bmu])
      for (i in 1:length(ind)) {
#        points(proto.vectors[ind[i],1], proto.vectors[ind[i],2],pch=16, col="green",cex=2/(H.distance[ind[i],bmu]+1))
         points(proto.vectors[ind[i],1], proto.vectors[ind[i],2],pch=16, col="white", cex=1.8)
         text(proto.vectors[ind[i],1], proto.vectors[ind[i],2],format(H.distance[ind[i],bmu]),col="blue",font=2)
        old.cb <- proto.vectors[ind[i],]
        proto.vectors[ind[i],] <- proto.vectors[ind[i],] - h.Dx[ind[i],]
        points(proto.vectors[ind[i],1], proto.vectors[ind[i],2], pch=16, col="red")
        lines(c(old.cb[1], proto.vectors[ind[i],1]), c(old.cb[2], proto.vectors[ind[i],2]), col="blue")
        readline("Press Enter...")
      }
      print("End of update")
         # The "expression" command will not allow mixing of expressions and
         # numbers so they are 'hard-wired'
    
      if (pass == 0) {
        connect.Map(5, 5,proto.vectors,  main=expression(paste("Input space - ",sigma^2,
            " = 17.64   ",alpha, " = 0.5")))
      } else if (pass == 1) {
        connect.Map(5, 5, proto.vectors,  main=expression(paste("Input space - ",sigma^2,
            " = 11.56   ",alpha, " = 0.5")))
      } else if (pass == 2) {
        connect.Map(5, 5, proto.vectors,  main=expression(paste("Input space - ",sigma^2,
            " = 6.76   ",alpha, " = 0.5")))
      } else if (pass == 3) {
        connect.Map(5, 5,proto.vectors,  main=expression(paste("Input space - ",sigma^2,
            " = 3.24   ",alpha, " = 0.5")))
      } else if (pass == 4) {
        connect.Map(5, 5, proto.vectors,  main=expression(paste("Input space - ",sigma^2,
            " = 1.00   ",alpha, " = 0.5")))
      }
      points(data[,1],data[,2],col="red",pch="+")
   }
   return(proto.vectors)
}