1 | #region License Information
|
---|
2 | /* HeuristicLab
|
---|
3 | * Copyright (C) 2002-2011 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 Kruskal-Shepard 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 | } |
---|