Changeset 10112
 Timestamp:
 11/07/13 00:21:52 (11 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

branches/HeuristicLab.Analysis.AlgorithmBehavior/HeuristicLab.Analysis.AlgorithmBehavior.Analyzers/3.3/DistanceMatrixToPoints.cs
r10109 r10112 27 27 public static class DistanceMatrixToPoints { 28 28 /* 29 * Calculates a matrix of ndimensional points from the distance matrix dm as described in 29 * Calculates a matrix of ndimensional points from the distance matrix dm as described in 30 30 * http://math.stackexchange.com/questions/156161/findingthecoordinatesofpointsfromdistancematrix/423898#423898 31 31 * and 32 32 * http://stackoverflow.com/questions/10963054/findingthecoordinatesofpointsfromdistancematrix/17177833#17177833 33 * 33 * 34 34 */ 35 35 public static double[][] ConvertDistanceMatrixToPoints(double[][] dm, int k = 2) { … … 76 76 //based on R's cmdscale 77 77 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]; 80 81 double[] q; //eigenvalues 81 double[,] v; //eigenvectors 82 double[,] v; //eigenvectors 83 82 84 83 85 double[][] x = SquareMatrix(dm); 84 86 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 85 133 ChangeSignAndHalve(x); 86 134 87 135 //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++) { 90 138 b[i, j] = x[i][j]; 91 139 } 92 140 } 93 141 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); 95 143 if (!res) throw new Exception("Eigenvalue computation did not converge!"); 96 144 97 145 //TODO: this should also work without allocating memory for ev and evec 98 146 double[] ev = new double[k]; 99 double[][] evec = new double[ dm.Length][];147 double[][] evec = new double[n][]; 100 148 AllocArray(evec, k); 101 149 Array.Copy(q, q.Length  k, ev, 0, k); 102 150 for (int i = 0; i < k; i++) { 103 for (int j = 0; j < dm.Length; j++) {151 for (int j = 0; j < n; j++) { 104 152 evec[j][i] = v[j, i + (q.Length  k)]; 105 153 } … … 114 162 AllocArray(points, k); 115 163 for (int i = 0; i < k; i++) { 116 for (int j = 0; j < dm.Length; j++) {164 for (int j = 0; j < n; j++) { 117 165 points[j][i] = Math.Sqrt(ev[i]) * evec[j][i]; 118 166 } … … 130 178 for (int j = 0; j < n; j++) { 131 179 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; 132 193 } 133 194 }
Note: See TracChangeset
for help on using the changeset viewer.