Changeset 10112 for branches/HeuristicLab.Analysis.AlgorithmBehavior/HeuristicLab.Analysis.AlgorithmBehavior.Analyzers

Timestamp:

11/07/13 00:21:52 (11 years ago)

Author:

ascheibe

Message:

#1886 added constant adding for metric mds

File:

• branches/HeuristicLab.Analysis.AlgorithmBehavior/HeuristicLab.Analysis.AlgorithmBehavior.Analyzers/3.3/DistanceMatrixToPoints.cs (modified) (4 diffs)

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

@@ -27 +27 @@
   public static class DistanceMatrixToPoints {
     /*
-     *  Calculates a matrix of n-dimensional points from the distance matrix dm as described in 
+     *  Calculates a matrix of n-dimensional points from the distance matrix dm as described in
      *  http://math.stackexchange.com/questions/156161/finding-the-coordinates-of-points-from-distance-matrix/423898#423898
      *  and
      *  http://stackoverflow.com/questions/10963054/finding-the-coordinates-of-points-from-distance-matrix/17177833#17177833
-     *   
+     *
      */
     public static double[][] ConvertDistanceMatrixToPoints(double[][] dm, int k = 2) {
@@ -76 +76 @@
     //based on R's cmdscale
     public static double[][] MetricMDS(double[][] dm, int k = 2, bool add = false) {
-      double[][] points = new double[dm.Length][];
-      double[,] b = new double[dm.Length, dm.Length];
+      int n = dm.Length;
+      double[][] points = new double[n][];
+      double[,] b = new double[n, n];
       double[] q; //eigenvalues
-      double[,] v; //eigenvectors 
+      double[,] v; //eigenvectors
+
 
       double[][] x = SquareMatrix(dm);
       CenterMatrix(x);
+
+      // solve additive constant problem
+      if (add) {
+        int[] i = Enumerable.Range(0, n).ToArray();
+        int[] i2 = Enumerable.Range(0, n).Select(y => y + n).ToArray();
+
+        double[,] Z = new double[n * 2, n * 2];
+
+        for (int j = 0; j < i.Length; j++) {
+          Z[i2[j], i[j]] = -1.0;
+        }
+
+        for (int j = 0; j < n; j++) {
+          for (int l = 0; l < n; l++) {
+            Z[i[j], i2[l]] = -1.0 * x[j][l];
+          }
+        }
+
+        double[][] centeredD = DoubleMatrix(dm);
+        CenterMatrix(centeredD);
+        for (int j = 0; j < n; j++) {
+          for (int l = 0; l < n; l++) {
+            Z[i2[j], i2[l]] = centeredD[j][l];
+          }
+        }
+
+        double[] wr;
+        double[] wi;
+        double[,] vl;
+        double[,] vr;
+        bool ret = alglib.rmatrixevd(Z, 2 * n, 0, out wr, out wi, out vl, out vr);
+        double c = wr.Max();
+
+        x = new double[n][];
+        AllocArray(x, n);
+
+        for (int j = 0; j < n; j++) {
+          for (int l = 0; l < n; l++) {
+            if (j != l) {
+              x[j][l] = Math.Pow(dm[j][l] + c, 2);
+            }
+          }
+        }
+        CenterMatrix(x);
+      }
+
       ChangeSignAndHalve(x);
 
       //TODO: optimize memory consumption
-      for (int i = 0; i < dm.Length; i++) {
-        for (int j = 0; j < dm.Length; j++) {
+      for (int i = 0; i < n; i++) {
+        for (int j = 0; j < n; j++) {
           b[i, j] = x[i][j];
         }
       }
 
-      bool res = alglib.smatrixevd(b, dm.Length, 1, true, out q, out v);
+      bool res = alglib.smatrixevd(b, n, 1, true, out q, out v);
       if (!res) throw new Exception("Eigenvalue computation did not converge!");
 
       //TODO: this should also work without allocating memory for ev and evec
       double[] ev = new double[k];
-      double[][] evec = new double[dm.Length][];
+      double[][] evec = new double[n][];
       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++) {
+        for (int j = 0; j < n; j++) {
           evec[j][i] = v[j, i + (q.Length - k)];
         }
@@ -114 +162 @@
       AllocArray(points, k);
       for (int i = 0; i < k; i++) {
-        for (int j = 0; j < dm.Length; j++) {
+        for (int j = 0; j < n; j++) {
           points[j][i] = Math.Sqrt(ev[i]) * evec[j][i];
         }
@@ -130 +178 @@
         for (int j = 0; j < n; j++) {
           newA[i][j] = Math.Pow(a[i][j], 2.0);
+        }
+      }
+      return newA;
+    }
+
+    private static double[][] DoubleMatrix(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] = a[i][j] * 2.0;
         }
       }