batch.train.SOM <- function(trainlen, cb, Data, g.units, n.cases, Hex.distance, radius, n.var) {

   Dx2 <- t(2*Data)     # Take 2*distance to the datum for use in cosine law
      # A constant distance (length of each datum vector)
   Data.const <- drop(apply(Data*Data, 1, sum))
   for (t in 1:trainlen) {
      Dist <- matrix((apply(cb*cb, 1, sum)), g.units, n.cases) -  cb%*%Dx2
      min.dist <- apply(Dist, 2, min)
      min.index <- apply(Dist, 2, order)[1,]
           # Updating step: Compute the average of the data vactors
           # in the neighbourhood of a prototype (codebook) vector and use that
           # average as the new value for the prototype (codebook) vector.
           # The neighbouring prototype (codebook) vectors are also altered based
           # on their closeness to the first. The effect is controlled
           # by the neighborhood function.
           # This effect is is given by the H matrix.
           # So, for each prototype (codebook) vector the new weight vector is
           #
           #      m = sum_i (h_i * d_i) / sum_i (h_i), 
           # 
           # where i denotes the index of data vector.
           # The values of neighbourhood function h_i are the same
           # for all data vectors belonging to the Voronoi set of the
           # same prototype (codebook) vector, so we calculate a partition matrix P
           # with elements p_ij=1 if the BMU of data vector j is i.
      H <- exp(-Hex.distance/(2*radius[t]))         
      P <- matrix(0, g.units, n.cases)
      P[cbind(min.index, 1:n.cases)]<-1
         # The sum of vectors in each Voronoi set are calculated as
         # (P*Data) and the neighbourhood is taken into account by
         # calculating a weighted sum of the
         # Voronoi sum (H*). The "activation" matrix A is the
         # denominator of the equation above.
      S <- H%*%(P%*%Data)
      A <- H%*%matrix(rep(apply(P,1,sum), n.var), g.units, n.var)
      cb <- S*(1/A)
   }
   return(cb)
}