source: branches/FitnessLandscapeAnalysis/VRPProblemAnalyzer/Program.cs @ 7232

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.IO;

namespace VRPProblemAnalyzer {
  public static class Maths {
    public static double StandardDeviation(this IEnumerable<double> data) {
      double average = data.Average();
      var individualDeviations = data.Select(num => Math.Pow(num - average, 2.0));
      return Math.Sqrt(individualDeviations.Average());
    }
  }
  
  class Program {

    private static double GetDistance(double[,] vertices, int source, int target) {
      return Math.Sqrt(
        Math.Pow((vertices[source, 0] - vertices[target, 0]), 2) + 
        Math.Pow((vertices[source, 1] - vertices[target, 1]), 2));
    }

    private static double[] GetDistances(double[,] vertices) {
      List<double> result = new List<double>();

      int cities = vertices.Length / 2;
      for (int i = 0; i < cities; i++) {
        for (int j = 0; j < cities; j++) {
          if (i != j) {
            result.Add(GetDistance(vertices, i, j));
          }
        }
      }

      return result.ToArray();
    }

    private static double[] GetAverageDistances(double[,] vertices) {
      int cities = vertices.Length / 2;
      double[] result = new double[cities];
      
      for (int i = 0; i < cities; i++) {
        double distance = 0;
        for (int j = 0; j < cities; j++) {
          if (i != j) {
            distance += GetDistance(vertices, i, j);
          }
        }
        result[i] = distance / (cities - 1);
      }

      return result;
    }

    private static double GetAverageDistance(TSPLIBParser instance) {
      double result = 0;

      double[] distances = GetDistances(instance.Vertices);
      result = distances.Average();

      return result;
    }

    private static double GetDistanceHeterogenity(TSPLIBParser instance) {
      double result = 0;

      result = GetAverageDistances(instance.Vertices).StandardDeviation();

      return result;
    }

    private static double GetDistance(double[] cust1, double[] cust2) {
      return Math.Sqrt(
                 Math.Pow((cust1[0] - cust2[0]), 2) +
                 Math.Pow((cust1[1] - cust2[1]), 2));
    }

    private static double GetDistance(Cluster cluster1, Cluster cluster2) {
      List<double> distances = new List<double>();
      foreach (double[] cust1 in cluster1.Items) {
        foreach (double[] cust2 in cluster2.Items) {
          distances.Add(GetDistance(cust1, cust2));
        }
      }

      return distances.Min();

      //return GetDistance(cluster1.Center, cluster2.Center);
    }

    private static double GetDistance(List<Cluster> clusters, Cluster cluster) {
      List<double> distances = new List<double>();
      foreach (Cluster other in clusters) {
        if (other != cluster) {
          distances.Add(GetDistance(cluster, other));
        }
      }

      if (distances.Count > 0)
        return distances.Min();
      else
        return 0;
    }

    private static double GetCohesion(Cluster cluster) {
      double cohesion = 0;

      List<double> distances = new List<double>();
      foreach (double[] cust1 in cluster.Items) {
        //option 1
        foreach (double[] cust2 in cluster.Items) {
          if (cust1 != cust2) {
            distances.Add(GetDistance(cust1, cust2));            
          }
        }

        //option2
        //distances.Add(GetDistance(cluster.Center, cust1));
      }

      if (distances.Count > 0)
        cohesion = distances.Average();

      return cohesion;
    }

    private static double GetClustering(TSPLIBParser instance) {
      double result = 0;

      int cities = instance.Vertices.Length / 2 - 1;
      double[,] coordinates = new double[cities, 2];
      for (int i = 1; i <= cities; i++) {
        coordinates[i - 1, 0] = instance.Vertices[i, 0];
        coordinates[i - 1, 1] = instance.Vertices[i, 1];
      }

      double clusterFactor = 0;
      double couplingFactor = 1;
      double cohesionFactor = 1;

      double maxQuality = double.MinValue;

      int clusters = -1;
      for (int i = 1; i <= cities; i++) {
        List<Cluster> collection = KMeansClustering.Cluster(coordinates, i, 9);

        if (collection.Count > 0) {
          List<double> distances = new List<double>();
          List<double> cohesions = new List<double>();
          foreach (Cluster cluster in collection) {
            distances.Add(GetDistance(collection, cluster));
            cohesions.Add(GetCohesion(cluster));
          }

          double distance = distances.Average();
          double cohesion = cohesions.Average();

          double quality = clusterFactor * (distances.Count / (double)cities) + couplingFactor * distance - cohesionFactor * cohesion;
          if (quality > maxQuality) {
            maxQuality = quality;
            clusters = distances.Count;
          }
        }
      }

      if (clusters != 0)
        result = (cities / (double)clusters) / cities;

      return result;
    }

    private static double GetDemandSize(TSPLIBParser instance) {
      double result = 0;

      result = instance.Demands.Average();

      return result;
    }

    private static double GetDemandHeterogenity(TSPLIBParser instance) {
      double result = 0;

      result = instance.Demands.StandardDeviation();

      return result;
    }

    private static void Normalize(TSPLIBParser instance) {
      //normalize demands
      for (int i = 0; i < instance.Demands.Length; i++) {
          instance.Demands[i] = instance.Demands[i] / instance.Capacity;
      }

      //normalize coordinates
      //-find bounds
      double maxX = double.MinValue;
      double maxY = double.MinValue;

      int cities = instance.Vertices.Length / 2;
      for (int i = 0; i < cities; i++) {
        if (instance.Vertices[i, 0] > maxX)
          maxX = instance.Vertices[i, 0];

        if (instance.Vertices[i, 1] > maxY)
          maxY = instance.Vertices[i, 1];
      }

      //-normalize
      double factorX = 1;
      double factorY = 1;
      if (maxX > maxY) {
        factorX = 1 / maxX;
        factorY = (maxY / maxX) / maxY;
      } else {
        factorX = (maxX / maxY) / maxX;
        factorY = 1 / maxY;
      }

      for (int i = 0; i < cities; i++) {
        instance.Vertices[i, 0] *= factorX;
        instance.Vertices[i, 1] *= factorY;
      }
    }

    static void Main(string[] args) {
      if (args.Length < 1) {
        Console.Error.WriteLine("Not enough arguments, please specify directory with VRP instances.");
        Environment.Exit(-1);
      }
      var path = args[0];
      using (StreamWriter sw = new StreamWriter(Path.Combine(path, "analysis.csv"))) {
        sw.WriteLine("Instance;Customers;Clustering;DistanceAvg;DistanceStdev;DemandAvg;DemandStdev");
        
        string[] instances = Directory.GetFiles(path, "*.vrp", SearchOption.AllDirectories);
        foreach (string instance in instances) {
          TSPLIBParser parser = new TSPLIBParser(instance);
          parser.Parse();
          Normalize(parser);

          sw.Write(Path.GetFileNameWithoutExtension(instance));
          sw.Write(';');
          sw.Write(parser.Vertices.Length / 2 - 1);
          sw.Write(';');
          sw.Write(GetClustering(parser));
          sw.Write(';');
          sw.Write(GetAverageDistance(parser));
          sw.Write(';');
          sw.Write(GetDistanceHeterogenity(parser));
          sw.Write(';');
          sw.Write(GetDemandSize(parser));
          sw.Write(';');
          sw.Write(GetDemandHeterogenity(parser));
          sw.Write(';');
          sw.WriteLine();
        }
      }
    }
  }
}