Free cookie consent management tool by TermsFeed Policy Generator

Ignore:
Timestamp:
11/07/13 00:21:52 (11 years ago)
Author:
ascheibe
Message:

#1886 added constant adding for metric mds

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/HeuristicLab.Analysis.AlgorithmBehavior/HeuristicLab.Analysis.AlgorithmBehavior.Analyzers/3.3/DistanceMatrixToPoints.cs

    r10109 r10112  
    2727  public static class DistanceMatrixToPoints {
    2828    /*
    29      *  Calculates a matrix of n-dimensional points from the distance matrix dm as described in 
     29     *  Calculates a matrix of n-dimensional points from the distance matrix dm as described in
    3030     *  http://math.stackexchange.com/questions/156161/finding-the-coordinates-of-points-from-distance-matrix/423898#423898
    3131     *  and
    3232     *  http://stackoverflow.com/questions/10963054/finding-the-coordinates-of-points-from-distance-matrix/17177833#17177833
    33      *   
     33     *
    3434     */
    3535    public static double[][] ConvertDistanceMatrixToPoints(double[][] dm, int k = 2) {
     
    7676    //based on R's cmdscale
    7777    public static double[][] MetricMDS(double[][] dm, int k = 2, bool add = false) {
    78       double[][] points = new double[dm.Length][];
    79       double[,] b = new double[dm.Length, dm.Length];
     78      int n = dm.Length;
     79      double[][] points = new double[n][];
     80      double[,] b = new double[n, n];
    8081      double[] q; //eigenvalues
    81       double[,] v; //eigenvectors
     82      double[,] v; //eigenvectors
     83
    8284
    8385      double[][] x = SquareMatrix(dm);
    8486      CenterMatrix(x);
     87
     88      // solve additive constant problem
     89      if (add) {
     90        int[] i = Enumerable.Range(0, n).ToArray();
     91        int[] i2 = Enumerable.Range(0, n).Select(y => y + n).ToArray();
     92
     93        double[,] Z = new double[n * 2, n * 2];
     94
     95        for (int j = 0; j < i.Length; j++) {
     96          Z[i2[j], i[j]] = -1.0;
     97        }
     98
     99        for (int j = 0; j < n; j++) {
     100          for (int l = 0; l < n; l++) {
     101            Z[i[j], i2[l]] = -1.0 * x[j][l];
     102          }
     103        }
     104
     105        double[][] centeredD = DoubleMatrix(dm);
     106        CenterMatrix(centeredD);
     107        for (int j = 0; j < n; j++) {
     108          for (int l = 0; l < n; l++) {
     109            Z[i2[j], i2[l]] = centeredD[j][l];
     110          }
     111        }
     112
     113        double[] wr;
     114        double[] wi;
     115        double[,] vl;
     116        double[,] vr;
     117        bool ret = alglib.rmatrixevd(Z, 2 * n, 0, out wr, out wi, out vl, out vr);
     118        double c = wr.Max();
     119
     120        x = new double[n][];
     121        AllocArray(x, n);
     122
     123        for (int j = 0; j < n; j++) {
     124          for (int l = 0; l < n; l++) {
     125            if (j != l) {
     126              x[j][l] = Math.Pow(dm[j][l] + c, 2);
     127            }
     128          }
     129        }
     130        CenterMatrix(x);
     131      }
     132
    85133      ChangeSignAndHalve(x);
    86134
    87135      //TODO: optimize memory consumption
    88       for (int i = 0; i < dm.Length; i++) {
    89         for (int j = 0; j < dm.Length; j++) {
     136      for (int i = 0; i < n; i++) {
     137        for (int j = 0; j < n; j++) {
    90138          b[i, j] = x[i][j];
    91139        }
    92140      }
    93141
    94       bool res = alglib.smatrixevd(b, dm.Length, 1, true, out q, out v);
     142      bool res = alglib.smatrixevd(b, n, 1, true, out q, out v);
    95143      if (!res) throw new Exception("Eigenvalue computation did not converge!");
    96144
    97145      //TODO: this should also work without allocating memory for ev and evec
    98146      double[] ev = new double[k];
    99       double[][] evec = new double[dm.Length][];
     147      double[][] evec = new double[n][];
    100148      AllocArray(evec, k);
    101149      Array.Copy(q, q.Length - k, ev, 0, k);
    102150      for (int i = 0; i < k; i++) {
    103         for (int j = 0; j < dm.Length; j++) {
     151        for (int j = 0; j < n; j++) {
    104152          evec[j][i] = v[j, i + (q.Length - k)];
    105153        }
     
    114162      AllocArray(points, k);
    115163      for (int i = 0; i < k; i++) {
    116         for (int j = 0; j < dm.Length; j++) {
     164        for (int j = 0; j < n; j++) {
    117165          points[j][i] = Math.Sqrt(ev[i]) * evec[j][i];
    118166        }
     
    130178        for (int j = 0; j < n; j++) {
    131179          newA[i][j] = Math.Pow(a[i][j], 2.0);
     180        }
     181      }
     182      return newA;
     183    }
     184
     185    private static double[][] DoubleMatrix(double[][] a) {
     186      int n = a.Length;
     187      double[][] newA = new double[a.Length][];
     188
     189      for (int i = 0; i < n; i++) {
     190        newA[i] = new double[a.Length];
     191        for (int j = 0; j < n; j++) {
     192          newA[i][j] = a[i][j] * 2.0;
    132193        }
    133194      }
Note: See TracChangeset for help on using the changeset viewer.