Changeset 10108

Timestamp:

11/05/13 23:04:16 (11 years ago)

Author:

ascheibe

Message:

#1886 fixed metric MDS and updated unit test

Location:

branches/HeuristicLab.Analysis.AlgorithmBehavior

Files:

• AlgorithmBehaviorUnitTests/AlgorithmBehaviorUnitTests.csproj (modified) (1 diff)
• AlgorithmBehaviorUnitTests/DistanceMatrixToPointsTest.cs (modified) (4 diffs)
• HeuristicLab.Analysis.AlgorithmBehavior.Analyzers/3.3/DistanceMatrixToPoints.cs (modified) (1 diff)

branches/HeuristicLab.Analysis.AlgorithmBehavior/AlgorithmBehaviorUnitTests/AlgorithmBehaviorUnitTests.csproj

@@ -36 +36 @@
   </PropertyGroup>
   <ItemGroup>
+    <Reference Include="ALGLIB-3.7.0, Version=3.7.0.0, Culture=neutral, PublicKeyToken=ba48961d6f65dcec, processorArchitecture=MSIL">
+      <SpecificVersion>False</SpecificVersion>
+      <HintPath>..\..\..\trunk\sources\bin\ALGLIB-3.7.0.dll</HintPath>
+    </Reference>
     <Reference Include="HeuristicLab.Common-3.3">
       <HintPath>..\..\..\trunk\sources\bin\HeuristicLab.Common-3.3.dll</HintPath>

branches/HeuristicLab.Analysis.AlgorithmBehavior/AlgorithmBehaviorUnitTests/DistanceMatrixToPointsTest.cs

@@ -22 +22 @@
 using System;
 using HeuristicLab.Analysis.AlgorithmBehavior.Analyzers;
+using HeuristicLab.Common;
 using Microsoft.VisualStudio.TestTools.UnitTesting;
 
@@ -28 +29 @@
   public class DistanceMatrixToPointsTest {
     [TestMethod]
-    public void TestDistanceMatrixConversion() {
-      int nrOfPoints = 4;
+    public void TestDistanceMatrixConversionStatic() {
+      int nrOfPoints = 3;
       int dim = 2;
       double[][] orgPoints = new double[nrOfPoints][];
@@ -36 +37 @@
       double[][] newPoints = null;
 
-
       AllocArray(orgPoints, dim);
       AllocArray(orgDm, nrOfPoints);
       AllocArray(newDm, nrOfPoints);
-      SamplePoints(orgPoints);
+      //SamplePoints(orgPoints);
+      StaticPoints(orgPoints);
       CalculateDistanceMatrix(orgDm, orgPoints);
 
-      newPoints = DistanceMatrixToPoints.ConvertDistanceMatrixToPoints(orgDm);
+      newPoints = DistanceMatrixToPoints.MetricMDS(orgDm);
 
       CalculateDistanceMatrix(newDm, newPoints);
-      //TODO
+      Console.WriteLine("orgDm:");
+      PrintDM(orgDm);
+      Console.WriteLine("newDm:");
+      PrintDM(newDm);
 
+      for (int i = 0; i < orgDm.Length; i++) {
+        for (int j = 0; j < orgDm.Length; j++) {
+          double diff = orgDm[i][j] - newDm[i][j];
+          Assert.IsTrue(diff.IsAlmost(0.0));
+        }
+      }
+    }
+
+    private static void PrintDM(double[][] dm) {
+      for (int i = 0; i < dm.Length; i++) {
+        for (int j = 0; j < dm.Length; j++) {
+          Console.Write(dm[i][j] + " ");
+        }
+        Console.WriteLine();
+      }
     }
 
@@ -60 +79 @@
     }
 
+    private static void StaticPoints(double[][] points) {
+      points[0][0] = 2;
+      points[0][1] = 1;
+
+      points[1][0] = 5;
+      points[1][1] = 5;
+
+      points[2][0] = 3;
+      points[2][1] = 7;
+    }
+
     private static void CalculateDistanceMatrix(double[][] dm, double[][] points) {
       for (int i = 0; i < points.Length; i++) {
         for (int j = 0; j < points.Length; j++) {
-          dm[i][j] = CalcEuclideanDistance(points[i], points[j]);
+          dm[i][j] = points[i].EuclideanDistance(points[j]);
         }
       }
-    }
-
-    private static double CalcEuclideanDistance(double[] p1, double[] p2) {
-      double result = 0.0;
-      for (int i = 0; i < p1.Length; i++) {
-        result += Math.Pow(p1[i] - p2[i], 2);
-      }
-      return Math.Sqrt(result);
     }
 

branches/HeuristicLab.Analysis.AlgorithmBehavior/HeuristicLab.Analysis.AlgorithmBehavior.Analyzers/3.3/DistanceMatrixToPoints.cs

@@ -72 +72 @@
     }
 
-    public static double[][] MetricMDS(double[][] dm) {
+    //based on R's cmdscale
+    public static double[][] MetricMDS(double[][] dm, int k = 2, bool add = false) {
       double[][] points = new double[dm.Length][];
-      double[,] a = new double[dm.Length, dm.Length];
       double[,] b = new double[dm.Length, dm.Length];
-      double[] ai = new double[dm.Length];
-      double[] aj = new double[dm.Length];
-      double aavg = 0.0;
-      double[] q = new double[dm.Length];
-      double[,] v = new double[dm.Length, dm.Length];
+      double[] q; //eigenvalues
+      double[,] v; //eigenvectors 
 
+      double[][] x = SquareMatrix(dm);
+      CenterMatrix(x);
+      ChangeSignAndHalve(x);
 
+      //TODO: optimize memory consumption
       for (int i = 0; i < dm.Length; i++) {
         for (int j = 0; j < dm.Length; j++) {
-          a[i, j] = -0.5 * Math.Pow(dm[i][j], 2);
-          aavg += a[i, j];
+          b[i, j] = x[i][j];
         }
       }
 
-      aavg /= dm.Length * dm.Length;
+      bool res = alglib.smatrixevd(b, dm.Length, 1, true, out q, out v);
+      if (!res) throw new Exception("Eigenvalue computation did not converge!");
 
-      for (int i = 0; i < dm.Length; i++) {
-        ai[i] = dm[i].Average();
-      }
-
-      for (int i = 0; i < dm.Length; i++) {
+      //TODO: this should also work without allocating memory for ev and evec
+      double[] ev = new double[k];
+      double[][] evec = new double[dm.Length][];
+      AllocArray(evec, k);
+      Array.Copy(q, q.Length - k, ev, 0, k);
+      for (int i = 0; i < k; i++) {
         for (int j = 0; j < dm.Length; j++) {
-          aj[i] += dm[j][i];
-        }
-        aj[i] /= dm.Length;
-      }
-
-      for (int i = 0; i < dm.Length; i++) {
-        for (int j = 0; j < dm.Length; j++) {
-          b[i, j] = a[i, j] - ai[i] - aj[j] + aavg;
+          evec[j][i] = v[j, i + (q.Length - k)];
         }
       }
 
-      //TODO: check b for negative values and make it positiv
-
-      double[] tau;
-      double[,] r;
-      alglib.rmatrixqr(ref b, dm.Length, dm.Length, out tau);
-      alglib.rmatrixqrunpackr(b, dm.Length, dm.Length, out r);
-
-      bool res = alglib.smatrixevd(r, dm.Length, 1, true, out q, out v);
-      if (!res) throw new Exception("Eigenvalue computation did not converge!");
-
-      int zeroCnt = q.Count(x => x.IsAlmost(0) || x < 0.0);
-      for (int i = 0; i < dm.Length; i++) {
-        points[i] = new double[dm.Length - zeroCnt];
+      double k1 = SumIfLZero(ev);
+      if (k1 < k) {
+        throw new Exception("Zero-eigenvalues detected. This leads to a degenerate point set. Use constants. ");
+        //TODO: handling of this case; implement adding of constants
       }
 
-      int pi = 0;
-      for (int i = 0; i < dm.Length; i++) {
-        if (!q[i].IsAlmost(0.0) && q[i] > 0.0) {
-          for (int j = 0; j < dm.Length; j++) {
-            points[j][pi] = Math.Sqrt(q[i]) * v[j, i];
-          }
-          pi++;
+      AllocArray(points, k);
+      for (int i = 0; i < k; i++) {
+        for (int j = 0; j < dm.Length; j++) {
+          points[j][i] = Math.Sqrt(ev[i]) * evec[j][i];
         }
       }
+      return points;
+    }
 
-      return points;
+    //TODO: refactor the following methods into something sane
+    private static double[][] SquareMatrix(double[][] a) {
+      int n = a.Length;
+      double[][] newA = new double[a.Length][];
+
+      for (int i = 0; i < n; i++) {
+        newA[i] = new double[a.Length];
+        for (int j = 0; j < n; j++) {
+          newA[i][j] = Math.Pow(a[i][j], 2.0);
+        }
+      }
+      return newA;
+    }
+
+    //based on R's DoubleCentre
+    private static void CenterMatrix(double[][] a) {
+      int n = a.Length;
+
+      //reduce lines by line avg
+      for (int i = 0; i < n; i++) {
+        double sum = 0;
+        for (int j = 0; j < n; j++) sum += a[i][j];
+        sum /= n;
+        for (int j = 0; j < n; j++) a[i][j] -= sum;
+      }
+
+      //reduce cols by col avg
+      for (int j = 0; j < n; j++) {
+        double sum = 0;
+        for (int i = 0; i < n; i++) sum += a[i][j];
+        sum /= n;
+        for (int i = 0; i < n; i++) a[i][j] -= sum;
+      }
+    }
+
+    private static void ChangeSignAndHalve(double[][] a) {
+      int n = a.Length;
+
+      for (int i = 0; i < n; i++) {
+        for (int j = 0; j < n; j++) {
+          a[i][j] = (-1.0 * a[i][j]) / 2;
+        }
+      }
+    }
+
+    private static double SumIfLZero(double[] a) {
+      return a.Where(x => x > 0.0 && !x.IsAlmost(0.0)).Sum();
+    }
+
+    private static void AllocArray(double[][] arr, int size) {
+      for (int i = 0; i < arr.Length; i++) {
+        arr[i] = new double[size];
+      }
     }