-------------------------------------------------------------------------
This file gives pseudo code for the algorithm presented in:

  "Genetic algorithm with deterministic crossover for vector quantization"
   Pattern Recognition Letters, 2000, Vol.21, No.1, pp.61-68

Updated: 7 December, 1999
Pasi Franti
franti@cs.uef.fi

-------------------------------------------------------------------------
Preface:
-------------------------------------------------------------------------
The following pseudo code gives a sketch of the implementation.
It is not complete, not very exact, and several less relevant
lower level routines are left out. The pseudo code aims at
readability and therefore few less important optimizations are
left out. The real source code uses several separate modules and
lots of experimental stuff, and therefore it is not available
in public. Questions and feedback by email are always welcome.


-------------------------------------------------------------------------
Data structures:
-------------------------------------------------------------------------

X[N] = Training set of N vectors.
P[N] = Partition of the training set, described as N pointers from X to C.
C[k] = Codebook of M vectors.
         C[j].vector = Code vector
         C[j].size   = Cluster size

S = Size of the population
T = Number of iterations
Solution refers to the pair (P,C).


-------------------------------------------------------------------------
Algorithm:
-------------------------------------------------------------------------


GeneticAlgorithm(X)
{
  GenerateInitialSolutions(X,S);
  SortSolutions;
  REPEAT T times
     {
     GenerateNewSolutions;
     SortSolutions;
     StoreBestSolution;
     }
  OutputBestCodebook;
}


------------------------- Initializations ---------------------------


GenerateInitialSolutions(X,S)
{
  FOR i=1 TO S DO
    {
    C_i <- GenerateRandomCodebook(X);
    P_i <- GenerateOptimalPartition(X,C_i);
    }
}


GenerateRandomCodebook(X)
{
  FOR j=1 TO k DO
    {
    /* The vector can be accepted as such but in our algorithm	*/
    /* we check that the same vector won't be chosen twice.     */
    i <- SelectRandomVector(X);
    C[j].vector <- X[i];
    C[j].size <- 0;
    }
  return C;
}


GenerateOptimalPartition(X,C)
{
  FOR i=1 TO N DO
     {
     j <- FindNearestVector(X[i],C);
     P[i] <- j;
     C[j].size <- C[j].size + 1;
     }
  return P;
}


------------------------- Main routines for GA ---------------------------


GenerateNewSolutions
{
  CalculateSurvivors(S);
  a <- 1; b <- 1; s <- 1;
  REPEAT S times
     {
     (a,b) <- SelectNextPair(a,b);
     (Cnew_s,Pnew_s) <- CrossSolutions(C_a, P_a, C_b, P_b);
     MutateSolution(Cnew_s);
     IterateByKMeans(Cnew_s, Pnew_s);
     s <- s + 1;
     }
  return {Cnew,Pnew};
}


CalculateSurvivors(S)
{
  s <- 1;
  WHILE s*(s-1)/2 < S DO
    {
    s <- s + 1;
    }
  NumberOfSurvivors = 1;
  CrossSetSize = s;
}


SelectNextPair(a,b)
/* Here a and b describe the previous pair */
{
  b <- b + 1;
  IF b=CrossSetSize THEN
     {
     a <- a + 1;
     b <- a + 1;
     }
  return (a,b);
}


CrossSolutions(C_1, P_1, C_2, P_2)
{
  Cnew <- CombineCentroids(C_1, C_2);
  Pnew <- CombinePartitions(P_1, P_2, C_1, C_2);
  Cnew <- UpdateCentroids(Pnew);
  RemoveEmptyClusters(Cnew, Pnew);
  PerformPNN(Cnew, Pnew);
  return (Cnew, Pnew);
}


CombineCentroids(C_1, C_2)
{
  i <- 1;
  FOR j=1 TO BookSize(C_1) DO
    {
    Cnew[i] <- C_1[j];
    i <- i + 1;
    }
  FOR j=1 TO BookSize(C_2) DO
    {
    Cnew[i] <- C_2[j];
    i <- i + 1;
    }
  return Cnew;
}


CombinePartitions(P_1, P_2, C_1, C_2)
{
  FOR i=1 TO N DO
     {
     p1 <- P_1[i];
     p2 <- P_2[i];
     d1 <- Distance(X[i], C_1[p1]);
     d2 <- Distance(X[i], C_2[p2]);
     IF d1 < d2 THEN Pnew[i] <- p1;
                ELSE Pnew[i] <- p2 + BookSize(C_1);
     /* The indices from P_1 will be in range [1..BookSize] but the  */
     /* indices from P_2 will be in range [Booksize+1..2*Booksize].  */
     }
  return Pnew;
}


UpdateCentroids(P)
{
  /* First we calculate the sum vectors for each partition.   */
  /* Remember also to initialize the counters by zero values. */
  FOR i=1 TO N DO
    {
    j <- P[i];
    Sum[j] <- Sum[j] + X[i];
    Count[j] <- Count[j] + 1;
    }

  /* Then we calculate the centroids as the average vectors */
  FOR j=1 TO k DO
    {
    /* Remember to check that Count is non-zero */
    C[j].vector <- Sum[j] / Count[j];
    C[j].size <- Count[j];
    }
  return C;
}


RemoveEmptyClusters(C,P)
{
  FOR j=1 TO BookSize(C) DO
     {
     /* Empty cluster is removed by moving the last entry in the    */
     /* codebook to the place of the empty vector. The last place   */
     /* is then removed. The partition indices must also be updated */
     /* so that for all P[i]=last -> P[i]=j.                        */
     IF C[j].size=0 AND BookSize(C)>M THEN
        {
        last = BookSize(C);
        C[j] <- C[last];
        RemoveVector(C,last);
        JoinPartitions(P,j,last);
        }
     }
}


MutateSolution(C)
{
  i <- SelectRandomVector(X);
  j <- SelectRandomVector(C);
  C[j].vector <- X[i];
}


------------------------------ PNN routines ------------------------------


PerformPNN(C,P)
{
  Q[BookSize(C)]: Nearest Neighbor Pointers;
      Q[j].nearest  = Index of the nearest neighbor cluster.
      Q[j].distance = Cost of merging the two clusters.
      Q[j].recalculate = Flag indicating if the pointer must be updated.

  FOR j=1 TO BookSize(C) DO
     {
     Q[j] <- FindNearestNeighbor(C,j);
     }
  WHILE BookSize(C) > k DO
     {
     a <- FindMinimumDistance(C,Q);
     b <- Q[a].nearest;
     MergeVectors(C,P,Q,a,b);
     UpdatePointers(C,Q);
     }
}


FindNearestNeighbor(C,a)
{
  q: Nearest Neighbor Pointer structure;

  q.nearest <- 1;
  q.distance <- Infinite;
  q.recalculate <- NO;
  FOR j=1 TO BookSize(C) DO
     {
     d <- MergeDistortion(C[a],C[j]);
     IF a<>j AND d<q.distance THEN
        {
        q.nearest <- j;
        q.distance <- d;
        }
    }
  return q;
}


FindMinimumDistance(C,Q)
{
  MinDist <- Infinite;
  MinIndex <- 1;
  FOR j=1 TO BookSize(C) DO
    {
    IF Q[j].distance<MinDist THEN
       {
       MinIndex <- j;
       MinDist <- Q[j].distance;
       }
    }
  return MinIndex;
}


MergeVectors(C,P,Q,a,b)
{
  IF a>b THEN Swap(a,b); /* So that a is smaller index. */
  last = BookSize(C);
  MarkClustersForRecalculation(C,Q,a,b);
  C[a].vector <- Centroid(C[a].vector,C[b].vector);
  JoinPartitions(P,C,a,b);
  FillEmptyPosition(C,Q,b,last);
  DecreaseBookSize(C);
}


MarkClustersForRecalculation(C,Q,a,b)
{
  FOR j=1 TO j<BookSize(C) DO
     {
     IF Q[i].nearest=a OR Q[i].nearest=b THEN Q[i].Recalculate = YES;
                                         ELSE Q[i].Recalculate = NO;
     }
}


JoinPartitions(P,C,a,b)
/* The partitions a and b are joined so that all vectors will be */
/* in a and cluster b will be empty. Pointers are updated.       */
{
  FOR i=1 TO N DO
    {
    IF P[i]=b THEN P[i]<-a;
    }
  C[a].size <- C[a].size + C[b].size;
}


FillEmptyPosition(C,Q,b,last)
/* Merging two vectors create empty place in the codebook. In this      */
/* routine the empty place is filled by the last entry in the codebook. */
{
  IF b<>last THEN
     {
     C[b] <- C[last];
     Q[b] <- Q[last];
     FOR j=1 TO j<BookSize(C) DO
        {
        IF Q[j].nearest=last THEN
           {
           Q[j].nearest <- b;
           }
        }
     }


UpdatePointers(C,Q)
{
  FOR j=1 TO j<BookSize(C) DO
     {
     IF Q[j].recalculate=YES THEN
        {
        Q[j] <- FindNearestNeighbor(C,j);
        Q[j].recalculate <- NO;
        }
     }
}


MergeDistortion(c1,c2):
  Calculates the merge cost of the cluster using equation (3) in the paper.

Centroid(c1,c2):
  Calculates the weighted average of the two input vectors.

BookSize(C):
  Returns the number of vectors in the (current) code book C.

DecreaseBookSize(C):
  Decreases the counter (the size of the codebook) by one.


------------------------------ K-means routines ------------------------------


IterateByKMeans(C,P)
{
  Active[k]: Boolean flag describing the activity of code vector.
  ChangedList[k]: List of active code vectors.

  iterations <- 0;
  /* Note that you must put a dummy in the ChangedList to get going. */
  WHILE iterations<MaxNumberOfIterations AND NotEmpty(ChangeList) DO
     {
     CopyCodebook(C,C_prev);
     C <- CalculateCentroids(P);
     FillEmptyPartitions(X,C,P); /* Not really necessary with GA. */
     DetectChangedCodeVectors(C_prev,C,Active,ChangeList);
     P <- ReducedSearchPartition(X,C,P,Active,ChangeList);
     iterations <- iterations + 1;
     }
}


DetectChangedCodeVectors(C_prev,C_new,Active,ChangedList)
{
  MakeEmpty(ChangedList);
  FOR j=1 TO j<BookSize(C_prev) DO
     {
     Active[j] = NO;
     IF C_prev[j]<>C_new[j] THEN
        {
        AddList(ChangedList,j);
        Active[j] = YES;
        }
     }
}


ReducedSearchPartition(X,C,P,Active,ChangeList)
{
  FOR i=1 TO N DO
     {
     j <- P[i]; /* Current partition of the training vector. */
     IF Active[j] THEN /* Code vector changed: perform full search. */
        {
        P[i] <- FindNearestVector(X[i],C);
        }
     ELSE /* Code vector not changed: search from the active vectors only */
        {
        P[i] <- FindNearestInSet(X[i],C,ChangedList);
        }
     }
  return P;
}


FindNearestInSet(vector,C,ChangedList)
{
  MinDist <- Infinite;
  MinIndex <- 1;

  FOR i=1 TO Size(ChangedList) DO
     {
     j <- ChangedList[i]; /* Changed list is actually set of pointers. */
     d <- Distance(vector,C[j]);
     IF d<MinDist THEN
        {
        MinIndex <- j;
        MinDist <- d;
        }
     }
  return MinIndex;
}


CalculateCentroids(P):
  Calculates partition centroids and uses them as codevectors.

NotEmpty(L):
  Returns TRUE if there are any members in the list.

MakeEmpty(L):
  Makes the list empty.