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

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

namespace VRPProblemAnalyzer {
  
  public class Program {

    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(Utils.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 += Utils.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;
      }
      instance.Capacity = 1.0;

      //normalize coordinates
      //-find bounds
      var cities = Utils.MatrixToPointList(instance.Vertices);
      var minX = cities.Min(c => c.X);
      var minY = cities.Min(c => c.Y);
      var maxX = cities.Max(c => c.X);
      var maxY = cities.Max(c => c.Y);
      var rangeX = maxX - minX;
      var rangeY = maxY - minY;
      var factor = Math.Sqrt(0.5)/Math.Max(rangeX, rangeY);

      for (int i = 0; i < instance.Vertices.GetLength(0); i++) {
        instance.Vertices[i, 0] = (cities[i].X-minX)*factor;
        instance.Vertices[i, 1] = (cities[i].Y-minY)*factor;
      }
    }

    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;GeographicExcentricity;DistanceExcentricity;DistanceDemandExcenetricity");
        
        string[] instances = Directory.GetFiles(path, "*.vrp", SearchOption.AllDirectories);
        foreach (string instance in instances) {
          TSPLIBParser parser = new TSPLIBParser(instance);
          parser.Parse();
          Normalize(parser);

          string instanceName = Path.GetFileNameWithoutExtension(instance);

          string solutionName = Path.Combine(Path.GetDirectoryName(instance), Path.GetFileNameWithoutExtension(instance) + ".opt");
          SolutionParser solution = new SolutionParser(solutionName);
          solution.Parse();

          sw.Write(instanceName);
          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.Write(DepotExcentricityCalculator.Geographic(parser.Vertices));
          sw.Write(';');
          sw.Write(DepotExcentricityCalculator.DistanceCentroid(parser.Vertices));
          sw.Write(';');
          sw.Write(DepotExcentricityCalculator.DemandDistanceCentroid(parser.Vertices, parser.Demands));
          sw.Write(';');
          sw.WriteLine();


          Image picture = PictureGenerator.GeneratePicture(parser, solution);
          picture.Save(Path.Combine(path, instanceName + ".png"));
        }
      }
      //Console.WriteLine("Done. Press Enter...");
      //Console.ReadLine();
    }
  }
}