Changeset 5648 for branches/QAP/HeuristicLab.Problems.QuadraticAssignment.Views/3.3/MultidimensionalScaling.cs

Timestamp:

03/10/11 01:49:10 (13 years ago)

Author:

abeham

Message:

#1330

• Unified QAP visualization in solution and problem view
• Fixed bug in gradient-descent gradient calculation when performing multidimensional scaling
• Extended QAPLIB parsers to cover some file format variations
• Added unit tests to check if all QAPLIB instances import without error
• Changed BestKnownSolution to be an OptionalValueParameter

Location:

branches/QAP

Files:

• . (modified) (1 prop)
• HeuristicLab.Problems.QuadraticAssignment.Views/3.3/MultidimensionalScaling.cs (modified) (5 diffs)

branches/QAP

@@ -1 +1 @@
 *.suo
+TestResults

@@ -1 +1 @@
 *.suo
+TestResults

branches/QAP/HeuristicLab.Problems.QuadraticAssignment.Views/3.3/MultidimensionalScaling.cs

@@ -23 +23 @@
 using HeuristicLab.Data;
 
-namespace HeuristicLab.Problems.QuadraticAssignment {
+namespace HeuristicLab.Problems.QuadraticAssignment.Views {
   public static class MultidimensionalScaling {
-
-    public static DoubleMatrix Classic(DoubleMatrix distances, out double error) {
-      if (distances == null) throw new ArgumentNullException("distances");
-      if (distances.Rows != distances.Columns) throw new ArgumentException("Distance matrix must be a square matrix.", "distances");
-      error = 0.0;
-
-      int dimension = distances.Rows;
-      if (dimension == 1) return new DoubleMatrix(new double[,] { { 0, 0 } });
-      else if (dimension == 2) return new DoubleMatrix(new double[,] { { 0, distances[0, 1] } });
-
-      double[,] A = new double[dimension + 1, dimension + 1];
-      for (int i = 1; i < dimension; i++) {
-        for (int j = i + 1; j < dimension + 1; j++) {
-          A[i, j] = -0.5 * distances[i - 1, j - 1] * distances[i - 1, j - 1];
-          A[j, i] = A[i, j];
-          A[0, 0] += A[i, j] / (dimension * dimension);
-          double avg = A[i, j] / (double)dimension;
-          A[0, j] += avg;
-          A[0, i] += avg;
-          A[i, 0] += avg;
-          A[j, 0] += avg;
-        }
-      }
-      double[,] B = new double[dimension, dimension];
-      for (int i = 0; i < dimension - 1; i++) {
-        for (int j = i + 1; j < dimension; j++) {
-          B[i, j] = A[i + 1, j + 1] - A[0, j + 1] - A[i + 1, 0] + A[0, 0];
-        }
-      }
-      double[] eigenvalues;
-      double[,] eigenvectors;
-      alglib.smatrixevd(B, dimension, 1, true, out eigenvalues, out eigenvectors);
-
-      DoubleMatrix coordinates = new DoubleMatrix(dimension, 2);
-      for (int i = 0; i < dimension; i++) {
-        for (int j = 0; j < dimension; j++) {
-          coordinates[i, 0] += eigenvectors[eigenvectors.GetLength(0) - 1, j]
-                             * distances[i, j];
-          coordinates[i, 1] += eigenvectors[eigenvectors.GetLength(0) - 2, j]
-                             * distances[i, j];
-        }
-      }
-
-      return coordinates;
-    }
 
     /// <summary>
@@ -77 +32 @@
     /// </summary>
     /// <param name="distances">A symmetric NxN matrix that specifies the distances between each element i and j. Diagonal elements are ignored.</param>
-    /// <param name="error">Returns the error between the fit distances and the actual distances.</param>
+    /// <param name="stress">Returns the stress between the fit distances and the actual distances.</param>
     /// <returns>A Nx2 matrix where the first column represents the x- and the second column the y coordinates.</returns>
-    public static DoubleMatrix Metric(DoubleMatrix distances, out double error) {
+    public static DoubleMatrix MetricByDistance(DoubleMatrix distances, out double stress) {
       if (distances == null) throw new ArgumentNullException("distances");
       if (distances.Rows != distances.Columns) throw new ArgumentException("Distance matrix must be a square matrix.", "distances");
-      error = 0.0;
+      stress = 0.0;
 
       int dimension = distances.Rows;
@@ -94 +49 @@
         coordinates[i, 1] = 10 * Math.Sin(rad * i);
       }
-      double epsg = 0.0000000001;
+
+      double epsg = 1e-7;
       double epsf = 0;
       double epsx = 0;
-      int maxits = 1000;
+      int maxits = 0;
       alglib.mincgstate state = null;
       alglib.mincgreport rep;
 
-      for (int iterations = 0; iterations < 200; iterations++) {
+      for (int iterations = 0; iterations < 10; iterations++) {
         for (int i = 0; i < dimension; i++) {
           double[] c = new double[] { coordinates[i, 0], coordinates[i, 1] };
@@ -111 +67 @@
             alglib.mincgrestartfrom(state, c);
           }
-          alglib.mincgoptimize(state, func, null, new Info(coordinates, distances, i));
+          alglib.mincgoptimize(state, StressGradient, null, new Info(coordinates, distances, i));
           alglib.mincgresults(state, out c, out rep);
 
@@ -118 +74 @@
         }
       }
-      error = 0;
+      stress = CalculateNormalizedStress(dimension, distances, coordinates);
+      return coordinates;
+    }
+
+    private static void StressGradient(double[] x, ref double func, double[] grad, object obj) {
+      func = 0; grad[0] = 0; grad[1] = 0;
+      Info info = (obj as Info);
+      for (int i = 0; i < info.Coordinates.Rows; i++) {
+        double c = info.Distances[info.Row, i];
+        if (i != info.Row) {
+          double a = info.Coordinates[i, 0];
+          double b = info.Coordinates[i, 1];
+          func += Stress(x, c, a, b);
+          grad[0] += ((2 * x[0] - 2 * a) * Math.Sqrt(x[1] * x[1] - 2 * b * x[1] + x[0] * x[0] - 2 * a * x[0] + b * b + a * a) - 2 * c * x[0] + 2 * a * c) / Math.Sqrt(x[1] * x[1] - 2 * b * x[1] + x[0] * x[0] - 2 * a * x[0] + b * b + a * a);
+          grad[1] += ((2 * x[1] - 2 * b) * Math.Sqrt(x[1] * x[1] - 2 * b * x[1] + x[0] * x[0] - 2 * a * x[0] + b * b + a * a) - 2 * c * x[1] + 2 * b * c) / Math.Sqrt(x[1] * x[1] - 2 * b * x[1] + x[0] * x[0] - 2 * a * x[0] + b * b + a * a);
+        }
+      }
+    }
+
+    private static double Stress(double[] x, double distance, double xCoord, double yCoord) {
+      return Stress(x[0], x[1], distance, xCoord, yCoord);
+    }
+
+    private static double Stress(double x, double y, double distance, double otherX, double otherY) {
+      double d = Math.Sqrt((x - otherX) * (x - otherX)
+                         + (y - otherY) * (y - otherY));
+      return (d - distance) * (d - distance);
+    }
+
+    public static double CalculateNormalizedStress(int dimension, DoubleMatrix distances, DoubleMatrix coordinates) {
+      double stress = 0;
       for (int i = 0; i < dimension - 1; i++) {
         for (int j = i + 1; j < dimension; j++) {
           if (distances[i, j] != 0) {
-            double dx = coordinates[i, 0] - coordinates[j, 0];
-            double dy = coordinates[i, 1] - coordinates[j, 1];
-            double d = Math.Sqrt(dx * dx + dy * dy);
-            error += ((distances[i, j] - d) * (distances[i, j] - d)) / (distances[i, j] * distances[i, j]);
+            stress += Stress(coordinates[i, 0], coordinates[i, 1], distances[i, j], coordinates[j, 0], coordinates[j, 1])
+              / (distances[i, j] * distances[i, j]);
           }
         }
       }
-      return coordinates;
-    }
-
-    public static void func(double[] x, ref double func, double[] grad, object obj) {
-      func = 0; grad[0] = 0; grad[1] = 0;
-      Info info = (obj as Info);
-      for (int i = 0; i < x.Length; i++) {
-        if (i != info.Row) {
-          double d = (x[0] - info.Coordinates[i, 0])
-                   * (x[0] - info.Coordinates[i, 0])
-                   + (x[1] - info.Coordinates[i, 1])
-                   * (x[1] - info.Coordinates[i, 1]);
-          func += (d - info.Distances[info.Row, i]) * (d - info.Distances[info.Row, i]);
-          grad[0] += 2 * x[0] - 2 * info.Coordinates[i, 0];
-          grad[1] += 2 * x[1] - 2 * info.Coordinates[i, 1];
-        }
-      }
+      return stress;
     }
 

@@ -1 +1 @@
 *.suo
+TestResults