Free cookie consent management tool by TermsFeed Policy Generator

Ignore:
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:
3 edited

Legend:

Unmodified
Added
Removed
  • branches/HeuristicLab.Analysis.AlgorithmBehavior/AlgorithmBehaviorUnitTests/AlgorithmBehaviorUnitTests.csproj

    r10104 r10108  
    3636  </PropertyGroup>
    3737  <ItemGroup>
     38    <Reference Include="ALGLIB-3.7.0, Version=3.7.0.0, Culture=neutral, PublicKeyToken=ba48961d6f65dcec, processorArchitecture=MSIL">
     39      <SpecificVersion>False</SpecificVersion>
     40      <HintPath>..\..\..\trunk\sources\bin\ALGLIB-3.7.0.dll</HintPath>
     41    </Reference>
    3842    <Reference Include="HeuristicLab.Common-3.3">
    3943      <HintPath>..\..\..\trunk\sources\bin\HeuristicLab.Common-3.3.dll</HintPath>
  • branches/HeuristicLab.Analysis.AlgorithmBehavior/AlgorithmBehaviorUnitTests/DistanceMatrixToPointsTest.cs

    r10078 r10108  
    2222using System;
    2323using HeuristicLab.Analysis.AlgorithmBehavior.Analyzers;
     24using HeuristicLab.Common;
    2425using Microsoft.VisualStudio.TestTools.UnitTesting;
    2526
     
    2829  public class DistanceMatrixToPointsTest {
    2930    [TestMethod]
    30     public void TestDistanceMatrixConversion() {
    31       int nrOfPoints = 4;
     31    public void TestDistanceMatrixConversionStatic() {
     32      int nrOfPoints = 3;
    3233      int dim = 2;
    3334      double[][] orgPoints = new double[nrOfPoints][];
     
    3637      double[][] newPoints = null;
    3738
    38 
    3939      AllocArray(orgPoints, dim);
    4040      AllocArray(orgDm, nrOfPoints);
    4141      AllocArray(newDm, nrOfPoints);
    42       SamplePoints(orgPoints);
     42      //SamplePoints(orgPoints);
     43      StaticPoints(orgPoints);
    4344      CalculateDistanceMatrix(orgDm, orgPoints);
    4445
    45       newPoints = DistanceMatrixToPoints.ConvertDistanceMatrixToPoints(orgDm);
     46      newPoints = DistanceMatrixToPoints.MetricMDS(orgDm);
    4647
    4748      CalculateDistanceMatrix(newDm, newPoints);
    48       //TODO
     49      Console.WriteLine("orgDm:");
     50      PrintDM(orgDm);
     51      Console.WriteLine("newDm:");
     52      PrintDM(newDm);
    4953
     54      for (int i = 0; i < orgDm.Length; i++) {
     55        for (int j = 0; j < orgDm.Length; j++) {
     56          double diff = orgDm[i][j] - newDm[i][j];
     57          Assert.IsTrue(diff.IsAlmost(0.0));
     58        }
     59      }
     60    }
     61
     62    private static void PrintDM(double[][] dm) {
     63      for (int i = 0; i < dm.Length; i++) {
     64        for (int j = 0; j < dm.Length; j++) {
     65          Console.Write(dm[i][j] + " ");
     66        }
     67        Console.WriteLine();
     68      }
    5069    }
    5170
     
    6079    }
    6180
     81    private static void StaticPoints(double[][] points) {
     82      points[0][0] = 2;
     83      points[0][1] = 1;
     84
     85      points[1][0] = 5;
     86      points[1][1] = 5;
     87
     88      points[2][0] = 3;
     89      points[2][1] = 7;
     90    }
     91
    6292    private static void CalculateDistanceMatrix(double[][] dm, double[][] points) {
    6393      for (int i = 0; i < points.Length; i++) {
    6494        for (int j = 0; j < points.Length; j++) {
    65           dm[i][j] = CalcEuclideanDistance(points[i], points[j]);
     95          dm[i][j] = points[i].EuclideanDistance(points[j]);
    6696        }
    6797      }
    68     }
    69 
    70     private static double CalcEuclideanDistance(double[] p1, double[] p2) {
    71       double result = 0.0;
    72       for (int i = 0; i < p1.Length; i++) {
    73         result += Math.Pow(p1[i] - p2[i], 2);
    74       }
    75       return Math.Sqrt(result);
    7698    }
    7799
  • branches/HeuristicLab.Analysis.AlgorithmBehavior/HeuristicLab.Analysis.AlgorithmBehavior.Analyzers/3.3/DistanceMatrixToPoints.cs

    r10098 r10108  
    7272    }
    7373
    74     public static double[][] MetricMDS(double[][] dm) {
     74    //based on R's cmdscale
     75    public static double[][] MetricMDS(double[][] dm, int k = 2, bool add = false) {
    7576      double[][] points = new double[dm.Length][];
    76       double[,] a = new double[dm.Length, dm.Length];
    7777      double[,] b = new double[dm.Length, dm.Length];
    78       double[] ai = new double[dm.Length];
    79       double[] aj = new double[dm.Length];
    80       double aavg = 0.0;
    81       double[] q = new double[dm.Length];
    82       double[,] v = new double[dm.Length, dm.Length];
     78      double[] q; //eigenvalues
     79      double[,] v; //eigenvectors
    8380
     81      double[][] x = SquareMatrix(dm);
     82      CenterMatrix(x);
     83      ChangeSignAndHalve(x);
    8484
     85      //TODO: optimize memory consumption
    8586      for (int i = 0; i < dm.Length; i++) {
    8687        for (int j = 0; j < dm.Length; j++) {
    87           a[i, j] = -0.5 * Math.Pow(dm[i][j], 2);
    88           aavg += a[i, j];
     88          b[i, j] = x[i][j];
    8989        }
    9090      }
    9191
    92       aavg /= dm.Length * dm.Length;
     92      bool res = alglib.smatrixevd(b, dm.Length, 1, true, out q, out v);
     93      if (!res) throw new Exception("Eigenvalue computation did not converge!");
    9394
    94       for (int i = 0; i < dm.Length; i++) {
    95         ai[i] = dm[i].Average();
    96       }
    97 
    98       for (int i = 0; i < dm.Length; i++) {
     95      //TODO: this should also work without allocating memory for ev and evec
     96      double[] ev = new double[k];
     97      double[][] evec = new double[dm.Length][];
     98      AllocArray(evec, k);
     99      Array.Copy(q, q.Length - k, ev, 0, k);
     100      for (int i = 0; i < k; i++) {
    99101        for (int j = 0; j < dm.Length; j++) {
    100           aj[i] += dm[j][i];
    101         }
    102         aj[i] /= dm.Length;
    103       }
    104 
    105       for (int i = 0; i < dm.Length; i++) {
    106         for (int j = 0; j < dm.Length; j++) {
    107           b[i, j] = a[i, j] - ai[i] - aj[j] + aavg;
     102          evec[j][i] = v[j, i + (q.Length - k)];
    108103        }
    109104      }
    110105
    111       //TODO: check b for negative values and make it positiv
    112 
    113       double[] tau;
    114       double[,] r;
    115       alglib.rmatrixqr(ref b, dm.Length, dm.Length, out tau);
    116       alglib.rmatrixqrunpackr(b, dm.Length, dm.Length, out r);
    117 
    118       bool res = alglib.smatrixevd(r, dm.Length, 1, true, out q, out v);
    119       if (!res) throw new Exception("Eigenvalue computation did not converge!");
    120 
    121       int zeroCnt = q.Count(x => x.IsAlmost(0) || x < 0.0);
    122       for (int i = 0; i < dm.Length; i++) {
    123         points[i] = new double[dm.Length - zeroCnt];
     106      double k1 = SumIfLZero(ev);
     107      if (k1 < k) {
     108        throw new Exception("Zero-eigenvalues detected. This leads to a degenerate point set. Use constants. ");
     109        //TODO: handling of this case; implement adding of constants
    124110      }
    125111
    126       int pi = 0;
    127       for (int i = 0; i < dm.Length; i++) {
    128         if (!q[i].IsAlmost(0.0) && q[i] > 0.0) {
    129           for (int j = 0; j < dm.Length; j++) {
    130             points[j][pi] = Math.Sqrt(q[i]) * v[j, i];
    131           }
    132           pi++;
     112      AllocArray(points, k);
     113      for (int i = 0; i < k; i++) {
     114        for (int j = 0; j < dm.Length; j++) {
     115          points[j][i] = Math.Sqrt(ev[i]) * evec[j][i];
    133116        }
    134117      }
     118      return points;
     119    }
    135120
    136       return points;
     121    //TODO: refactor the following methods into something sane
     122    private static double[][] SquareMatrix(double[][] a) {
     123      int n = a.Length;
     124      double[][] newA = new double[a.Length][];
     125
     126      for (int i = 0; i < n; i++) {
     127        newA[i] = new double[a.Length];
     128        for (int j = 0; j < n; j++) {
     129          newA[i][j] = Math.Pow(a[i][j], 2.0);
     130        }
     131      }
     132      return newA;
     133    }
     134
     135    //based on R's DoubleCentre
     136    private static void CenterMatrix(double[][] a) {
     137      int n = a.Length;
     138
     139      //reduce lines by line avg
     140      for (int i = 0; i < n; i++) {
     141        double sum = 0;
     142        for (int j = 0; j < n; j++) sum += a[i][j];
     143        sum /= n;
     144        for (int j = 0; j < n; j++) a[i][j] -= sum;
     145      }
     146
     147      //reduce cols by col avg
     148      for (int j = 0; j < n; j++) {
     149        double sum = 0;
     150        for (int i = 0; i < n; i++) sum += a[i][j];
     151        sum /= n;
     152        for (int i = 0; i < n; i++) a[i][j] -= sum;
     153      }
     154    }
     155
     156    private static void ChangeSignAndHalve(double[][] a) {
     157      int n = a.Length;
     158
     159      for (int i = 0; i < n; i++) {
     160        for (int j = 0; j < n; j++) {
     161          a[i][j] = (-1.0 * a[i][j]) / 2;
     162        }
     163      }
     164    }
     165
     166    private static double SumIfLZero(double[] a) {
     167      return a.Where(x => x > 0.0 && !x.IsAlmost(0.0)).Sum();
     168    }
     169
     170    private static void AllocArray(double[][] arr, int size) {
     171      for (int i = 0; i < arr.Length; i++) {
     172        arr[i] = new double[size];
     173      }
    137174    }
    138175
Note: See TracChangeset for help on using the changeset viewer.