1  #region License Information


2  /* HeuristicLab


3  * Copyright (C) 20022011 Heuristic and Evolutionary Algorithms Laboratory (HEAL)


4  *


5  * This file is part of HeuristicLab.


6  *


7  * HeuristicLab is free software: you can redistribute it and/or modify


8  * it under the terms of the GNU General Public License as published by


9  * the Free Software Foundation, either version 3 of the License, or


10  * (at your option) any later version.


11  *


12  * HeuristicLab is distributed in the hope that it will be useful,


13  * but WITHOUT ANY WARRANTY; without even the implied warranty of


14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the


15  * GNU General Public License for more details.


16  *


17  * You should have received a copy of the GNU General Public License


18  * along with HeuristicLab. If not, see <http://www.gnu.org/licenses/>.


19  */


20  #endregion


21 


22  using System;


23  using HeuristicLab.Data;


24 


25  namespace HeuristicLab.Problems.QuadraticAssignment {


26  public static class MultidimensionalScaling {


27 


28  public static DoubleMatrix Classic(DoubleMatrix distances, out double error) {


29  if (distances == null) throw new ArgumentNullException("distances");


30  if (distances.Rows != distances.Columns) throw new ArgumentException("Distance matrix must be a square matrix.", "distances");


31  error = 0.0;


32 


33  int dimension = distances.Rows;


34  if (dimension == 1) return new DoubleMatrix(new double[,] { { 0, 0 } });


35  else if (dimension == 2) return new DoubleMatrix(new double[,] { { 0, distances[0, 1] } });


36 


37  double[,] A = new double[dimension + 1, dimension + 1];


38  for (int i = 1; i < dimension; i++) {


39  for (int j = i + 1; j < dimension + 1; j++) {


40  A[i, j] = 0.5 * distances[i  1, j  1] * distances[i  1, j  1];


41  A[j, i] = A[i, j];


42  A[0, 0] += A[i, j] / (dimension * dimension);


43  double avg = A[i, j] / (double)dimension;


44  A[0, j] += avg;


45  A[0, i] += avg;


46  A[i, 0] += avg;


47  A[j, 0] += avg;


48  }


49  }


50  double[,] B = new double[dimension, dimension];


51  for (int i = 0; i < dimension  1; i++) {


52  for (int j = i + 1; j < dimension; j++) {


53  B[i, j] = A[i + 1, j + 1]  A[0, j + 1]  A[i + 1, 0] + A[0, 0];


54  }


55  }


56  double[] eigenvalues;


57  double[,] eigenvectors;


58  alglib.smatrixevd(B, dimension, 1, true, out eigenvalues, out eigenvectors);


59 


60  DoubleMatrix coordinates = new DoubleMatrix(dimension, 2);


61  for (int i = 0; i < dimension; i++) {


62  for (int j = 0; j < dimension; j++) {


63  coordinates[i, 0] += eigenvectors[eigenvectors.GetLength(0)  1, j]


64  * distances[i, j];


65  coordinates[i, 1] += eigenvectors[eigenvectors.GetLength(0)  2, j]


66  * distances[i, j];


67  }


68  }


69 


70  return coordinates;


71  }


72 


73  /// <summary>


74  /// Performs the KruskalShepard algorithm and applies a gradient descent method


75  /// to fit the coordinates such that the difference between the fit distances


76  /// and the actual distances is minimal.


77  /// </summary>


78  /// <param name="distances">A symmetric NxN matrix that specifies the distances between each element i and j. Diagonal elements are ignored.</param>


79  /// <param name="error">Returns the error between the fit distances and the actual distances.</param>


80  /// <returns>A Nx2 matrix where the first column represents the x and the second column the y coordinates.</returns>


81  public static DoubleMatrix Metric(DoubleMatrix distances, out double error) {


82  if (distances == null) throw new ArgumentNullException("distances");


83  if (distances.Rows != distances.Columns) throw new ArgumentException("Distance matrix must be a square matrix.", "distances");


84  error = 0.0;


85 


86  int dimension = distances.Rows;


87  if (dimension == 1) return new DoubleMatrix(new double[,] { { 0, 0 } });


88  else if (dimension == 2) return new DoubleMatrix(new double[,] { { 0, distances[0, 1] } });


89 


90  DoubleMatrix coordinates = new DoubleMatrix(dimension, 2);


91  double rad = (2 * Math.PI) / coordinates.Rows;


92  for (int i = 0; i < dimension; i++) {


93  coordinates[i, 0] = 10 * Math.Cos(rad * i);


94  coordinates[i, 1] = 10 * Math.Sin(rad * i);


95  }


96  double epsg = 0.0000000001;


97  double epsf = 0;


98  double epsx = 0;


99  int maxits = 1000;


100  alglib.mincgstate state = null;


101  alglib.mincgreport rep;


102 


103  for (int iterations = 0; iterations < 200; iterations++) {


104  for (int i = 0; i < dimension; i++) {


105  double[] c = new double[] { coordinates[i, 0], coordinates[i, 1] };


106 


107  if (iterations == 0 && i == 0) {


108  alglib.mincgcreate(c, out state);


109  alglib.mincgsetcond(state, epsg, epsf, epsx, maxits);


110  } else {


111  alglib.mincgrestartfrom(state, c);


112  }


113  alglib.mincgoptimize(state, func, null, new Info(coordinates, distances, i));


114  alglib.mincgresults(state, out c, out rep);


115 


116  coordinates[i, 0] = c[0];


117  coordinates[i, 1] = c[1];


118  }


119  }


120  error = 0;


121  for (int i = 0; i < dimension  1; i++) {


122  for (int j = i + 1; j < dimension; j++) {


123  if (distances[i, j] != 0) {


124  double dx = coordinates[i, 0]  coordinates[j, 0];


125  double dy = coordinates[i, 1]  coordinates[j, 1];


126  double d = Math.Sqrt(dx * dx + dy * dy);


127  error += ((distances[i, j]  d) * (distances[i, j]  d)) / (distances[i, j] * distances[i, j]);


128  }


129  }


130  }


131  return coordinates;


132  }


133 


134  public static void func(double[] x, ref double func, double[] grad, object obj) {


135  func = 0; grad[0] = 0; grad[1] = 0;


136  Info info = (obj as Info);


137  for (int i = 0; i < x.Length; i++) {


138  if (i != info.Row) {


139  double d = (x[0]  info.Coordinates[i, 0])


140  * (x[0]  info.Coordinates[i, 0])


141  + (x[1]  info.Coordinates[i, 1])


142  * (x[1]  info.Coordinates[i, 1]);


143  func += (d  info.Distances[info.Row, i]) * (d  info.Distances[info.Row, i]);


144  grad[0] += 2 * x[0]  2 * info.Coordinates[i, 0];


145  grad[1] += 2 * x[1]  2 * info.Coordinates[i, 1];


146  }


147  }


148  }


149 


150  private class Info {


151  public DoubleMatrix Coordinates { get; set; }


152  public DoubleMatrix Distances { get; set; }


153  public int Row { get; set; }


154 


155  public Info(DoubleMatrix c, DoubleMatrix d, int r) {


156  Coordinates = c;


157  Distances = d;


158  Row = r;


159  }


160  }


161  }


162  } 
