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

#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) {
      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));
        }
      }

      //QR decomposition to get the upper part for smatrixevd
      double[] tau;
      double[,] r;
      alglib.rmatrixqr(ref m, dm.Length, dm.Length, out tau);
      alglib.rmatrixqrunpackr(m, 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];
      }

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

      return points;
    }

    public static double[][] MetricMDS(double[][] dm) {
      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];


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

      aavg /= dm.Length * dm.Length;

      for (int i = 0; i < dm.Length; i++) {
        ai[i] = dm[i].Average();
      }

      for (int i = 0; i < dm.Length; 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;
        }
      }

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

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

      return points;
    }
  }
}