source: branches/HeuristicLab.VRPEnhancements/HeuristicLab.Problems.VehicleRouting/3.4/Encodings/Potvin/Creators/Cluster.cs @ 14442

using System;
using System.Collections;
using System.Collections.Generic;
using HeuristicLab.Core;
using HeuristicLab.Problems.VehicleRouting.Interfaces;

namespace HeuristicLab.Problems.VehicleRouting.Encodings.Potvin {

  public class ClusterAlgorithm<C,CO> where C : Cluster<CO> where CO : ClusterElement {

    public static List<C> KMeans(IRandom random, List<CO> clusterElements, int k, double changeThreshold) {
      HashSet<int> initMeans = new HashSet<int>();
      int nextMean = -1;

      List<C> clusters = CreateCList();

      // (1) initialize each cluster with a random element as mean
      for (int i = 0; i < k && i < clusterElements.Count; i++) {
        C cluster = (C)Activator.CreateInstance(typeof(C));
        cluster.Id = i;

        do {
          nextMean = random.Next(0, clusterElements.Count);
        } while (initMeans.Contains(nextMean));
        initMeans.Add(nextMean);
        cluster.SetMean(clusterElements[nextMean]);
        clusters.Add(cluster);
      }

      // (2) repeat clustering until change rate is below threshold
      double changeRate = 0.0;

      do {
        int changes = KMeansRun(clusters, clusterElements);
        changeRate = (double)changes / clusterElements.Count;
      } while (changeRate > changeThreshold);


      // remove empty clusters
      clusters.RemoveAll(c => c.Elements.Count.Equals(0));

      return clusters;
    }

    private static int KMeansRun(List<C> clusters, List<CO> clusterElements) {
      int changes = 0;

      // clear clusters from previous assigned elements
      foreach (var c in clusters) {
        c.Elements.Clear();
      }

      // assign elements to currently most suitable clusters
      foreach (var e in clusterElements) {
        int optClusterIdx = 0;
        double optImpact = clusters[optClusterIdx].CalculateImpact(e);
        for (int i = 1; i < clusters.Count; i++) {
          double impact = clusters[i].CalculateImpact(e);
          if (impact < optImpact) { 
            optImpact = impact;
            optClusterIdx = i;
          }
        }
        if (clusters[optClusterIdx].AddObject(e)) {
          changes++;
        }
      }

      // update mean and variance
      foreach (var c in clusters) {
        c.CalculateMean();
        c.CalculateVariance();
      }

      return changes;
    }

    private static List<C> CreateCList() {
      var listType = typeof(List<>);
      var constructedListType = listType.MakeGenericType(typeof(C));
      return (List<C>)Activator.CreateInstance(constructedListType);
    }
  }

  #region Cluster
  public interface ICluster<T> where T : ClusterElement {
    void SetMean(T o);
    bool AddObject(T o);
    void CalculateMean();
    void CalculateVariance();
    double CalculateImpact(T e);
    double CalculateDistance(T e1, T e2);

  }
  public abstract class Cluster<T> : ICluster<T> where T : ClusterElement {
    public int Id;

    public List<T> Elements { get; private set; }

    public T Mean { get; set; }

    public double Variance { get; set; }

    protected Cluster() {
      Elements = new List<T>();
    }

    protected Cluster(int id) {
      Id = id;
      Elements = new List<T>();

    }

    public bool AddObject(T e) {
      Elements.Add(e);

      bool clusterChanged = e.ClusterId != Id;
      e.ClusterId = Id;
      return clusterChanged;
    }

    public void SetMean(T e) {
      Mean = e;
      Mean.ClusterId = Id;
    }

    public abstract void CalculateMean();

    public abstract double CalculateDistance(T e1, T e2);

    public virtual void CalculateVariance() {
      if (Mean == null)
        CalculateMean();

      Variance = 0.0;
      foreach (T e in Elements) {
        Variance += Math.Pow(CalculateDistance(Mean, e), 2);
      }
      Variance /= Elements.Count;
    }

    public virtual double CalculateImpact(T e) {
      if (Mean == null)
        CalculateMean();

      double newVariance = (Variance * Elements.Count + Math.Pow(CalculateDistance(Mean, e), 2)) / (Elements.Count + 1);
      return newVariance - Variance;
    }

  }
  public class SpatialDistanceCluster : Cluster<SpatialDistanceClusterElement> {
    public SpatialDistanceCluster() : base() {
    }

    public SpatialDistanceCluster(int id) : base(id) {
    }

    public override void CalculateMean() {
      int dimensions = (Mean != null) ? Mean.Coordinates.Length : (Elements.Count > 0) ? Elements[0].Coordinates.Length : 0;
      SpatialDistanceClusterElement mean = new SpatialDistanceClusterElement(-1, dimensions, Id);
      foreach (SpatialDistanceClusterElement e in Elements) {
        for (int i = 0; i < Mean.Coordinates.Length; i++) {
          mean.Coordinates[i] += (e.Coordinates[i] / Elements.Count);
        }
      }
      Mean = mean;
    }

    public override double CalculateDistance(SpatialDistanceClusterElement e1, SpatialDistanceClusterElement e2) {
      if (!e1.Coordinates.Length.Equals(e2.Coordinates.Length)) {
        throw new ArgumentException("Distance could not be calculated since number of dimensions is unequal.");
      }

      double distance = 0.0;
      for (int i = 0; i < e1.Coordinates.Length; i++) {
        distance += Math.Pow(e1.Coordinates[i] - e2.Coordinates[i], 2);
      }

      return Math.Sqrt(distance);
    }
  }

  public class TemporalDistanceCluster : Cluster<TemporalDistanceClusterElement> {
    public TemporalDistanceCluster() : base() {
    }

    public TemporalDistanceCluster(int id) : base(id) {
    }

    public override void CalculateMean() {
      TemporalDistanceClusterElement mean = new TemporalDistanceClusterElement(-1, Id);
      foreach (TemporalDistanceClusterElement e in Elements) {
        mean.ReadyTime += e.ReadyTime/Elements.Count;
        mean.DueTime += e.DueTime/Elements.Count;
      }
      Mean = mean;
    }

    public override double CalculateDistance(TemporalDistanceClusterElement e1, TemporalDistanceClusterElement e2) {
      double distance = 0.0;

      distance += Math.Pow(e1.ReadyTime - e2.ReadyTime, 2);
      distance += Math.Pow(e1.DueTime - e2.DueTime, 2);

      return Math.Sqrt(distance);
    }
  }

  #endregion Cluster

  #region ClusterElement
  public abstract class ClusterElement {
    public int Id { get; set; }
    public int ClusterId { get; set; }

    protected ClusterElement() {
      ClusterId = -1;
    }

    public ClusterElement(int id, int clusterId = -1) {
      Id = id;
      ClusterId = clusterId;
    }
  }
  public class SpatialDistanceClusterElement : ClusterElement {

    public double[] Coordinates { get; set; }

    public SpatialDistanceClusterElement() : base() {
    }

    public SpatialDistanceClusterElement(int id, int clusterId = -1) : base(id, clusterId) {
    }

    public SpatialDistanceClusterElement(int id, double[] coordinates, int clusterId = -1) : base(id, clusterId) {
      Coordinates = coordinates;
    }

    public SpatialDistanceClusterElement(int id, int dimensions, int clusterId = -1) : base(id, clusterId) { 
      Coordinates = new double[dimensions];
    }

  }

  public class TemporalDistanceClusterElement : ClusterElement {
    public double ReadyTime { get; set; }
    public double DueTime { get; set; }

    public TemporalDistanceClusterElement() : base() {
    }

    public TemporalDistanceClusterElement(int id, int clusterId = -1) : base(id, clusterId) {
    }

    public TemporalDistanceClusterElement(int id, double readyTime, double dueTime, int clusterId = -1) : base(id, clusterId) {
      ReadyTime = readyTime;
      DueTime = dueTime;
    }
  }

  #endregion ClusterElement
}