source: branches/HeuristicLab.Analysis.AlgorithmBehavior/HeuristicLab.Analysis.AlgorithmBehavior.Analyzers/3.3/DistanceMatrixToPoints.cs @ 10635

#region License Information
/* HeuristicLab
 * Copyright (C) 2002-2013 Heuristic and Evolutionary Algorithms Laboratory (HEAL)
 *
 * This file is part of HeuristicLab.
 *
 * HeuristicLab is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * HeuristicLab is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with HeuristicLab. If not, see <http://www.gnu.org/licenses/>.
 */
#endregion

using System;
using System.Linq;
using HeuristicLab.Common;

namespace HeuristicLab.Analysis.AlgorithmBehavior.Analyzers {
  public static class DistanceMatrixToPoints {
    /*
     *  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) {
      double[][] points = new double[dm.Length][];
      double[,] m = new double[dm.Length, dm.Length];
      double[] q = new double[dm.Length]; //eigenvalues
      double[,] v = new double[dm.Length, dm.Length]; //eigenvectors

      for (int i = 0; i < dm.Length; i++) {
        for (int j = 0; j < dm.Length; j++) {
          m[i, j] = 0.5 * (Math.Pow(dm[0][j], 2) + Math.Pow(dm[i][0], 2) - Math.Pow(dm[i][j], 2));
        }
      }

      bool res = alglib.smatrixevd(m, dm.Length, 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][];
      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++) {
          evec[j][i] = v[j, i + (q.Length - k)];
        }
      }

      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
      }

      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;
    }

    //based on R's cmdscale
    public static double[][] MetricMDS(double[][] dm, int k = 2, bool add = false) {
      int n = dm.Length;
      double[][] points = new double[n][];
      double[,] b = new double[n, n];
      double[] q; //eigenvalues
      double[,] v; //eigenvectors

      if (n < k)
        throw new ArgumentException("Distance matrix length must be greater than dimension", "k");

      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 < n; i++) {
        for (int j = 0; j < n; j++) {
          b[i, j] = x[i][j];
        }
      }

      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[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 < n; j++) {
          evec[j][i] = v[j, i + (q.Length - k)];
        }
      }

      int k1 = SumIfLZero(ev);
      if (k1 < k) {
        throw new Exception("Zero-eigenvalues detected. This leads to a degenerate point set. Use constants. ");
      }

      AllocArray(points, k);
      for (int i = 0; i < k; i++) {
        for (int j = 0; j < n; j++) {
          points[j][i] = Math.Sqrt(ev[i]) * evec[j][i];
        }
      }
      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;
    }

    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;
        }
      }
      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 int SumIfLZero(double[] a) {
      return a.Count(x => x > 0.0 && !x.IsAlmost(0.0));
    }

    private static void AllocArray(double[][] arr, int size) {
      for (int i = 0; i < arr.Length; i++) {
        arr[i] = new double[size];
      }
    }

    public static double[][] TransformToDistances(double[][] similarityMatrix) {
      double[][] dm = new double[similarityMatrix.Length][];

      for (int i = 0; i < dm.Length; i++) {
        dm[i] = new double[similarityMatrix.Length];
        for (int j = 0; j < dm.Length; j++) {
          dm[i][j] = Math.Sqrt(similarityMatrix[i][i] + similarityMatrix[j][j] - 2 * similarityMatrix[i][j]);
        }
      }

      return dm;
    }
  }
}