 Timestamp:
 03/27/17 17:27:03 (6 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/TSNEStatic.cs
r14785 r14788 57 57 using System.Collections.Generic; 58 58 using System.Linq; 59 using HeuristicLab.Analysis;60 59 using HeuristicLab.Collections; 61 60 using HeuristicLab.Common; 62 61 using HeuristicLab.Core; 63 using HeuristicLab.Data;64 using HeuristicLab.Optimization;65 62 using HeuristicLab.Persistence.Default.CompositeSerializers.Storable; 66 63 using HeuristicLab.Random; … … 68 65 namespace HeuristicLab.Algorithms.DataAnalysis { 69 66 [StorableClass] 70 public class TSNE<T> : DeepCloneable /*where T : class, IDeepCloneable*/ { 71 72 private const string IterationResultName = "Iteration"; 73 private const string ErrorResultName = "Error"; 74 private const string ErrorPlotResultName = "ErrorPlot"; 75 private const string ScatterPlotResultName = "Scatterplot"; 76 private const string DataResultName = "Projected Data"; 77 78 #region Properties 79 [Storable] 80 private IDistance<T> distance; 81 [Storable] 82 private int maxIter; 83 [Storable] 84 private int stopLyingIter; 85 [Storable] 86 private int momSwitchIter; 87 [Storable] 88 double momentum; 89 [Storable] 90 private double finalMomentum; 91 [Storable] 92 private double eta; 93 [Storable] 94 private IRandom random; 95 [Storable] 96 private ResultCollection results; 97 [Storable] 98 private Dictionary<string, List<int>> dataRowLookup; 99 [Storable] 100 private Dictionary<string, ScatterPlotDataRow> dataRows; 101 #endregion 102 103 #region Stopping 104 public volatile bool Running; // TODO 105 #endregion 106 107 #region HLConstructors & Cloning 108 [StorableConstructor] 109 protected TSNE(bool deserializing) { } 110 protected TSNE(TSNE<T> original, Cloner cloner) : base(original, cloner) { 111 distance = cloner.Clone(original.distance); 112 maxIter = original.maxIter; 113 stopLyingIter = original.stopLyingIter; 114 momSwitchIter = original.momSwitchIter; 115 momentum = original.momentum; 116 finalMomentum = original.finalMomentum; 117 eta = original.eta; 118 random = cloner.Clone(random); 119 results = cloner.Clone(results); 120 dataRowLookup = original.dataRowLookup.ToDictionary(entry => entry.Key, entry => entry.Value.Select(x => x).ToList()); 121 dataRows = original.dataRows.ToDictionary(entry => entry.Key, entry => cloner.Clone(entry.Value)); 67 public class TSNE<T> { 68 69 [StorableClass] 70 public sealed class TSNEState : DeepCloneable { 71 // initialized once 72 public IDistance<T> distance; 73 public IRandom random; 74 public double perplexity; 75 public bool exact; 76 public int noDatapoints; 77 public double finalMomentum; 78 public int momSwitchIter; 79 public int stopLyingIter; 80 public double theta; 81 public double eta; 82 public int newDimensions; 83 84 // for approximate version: sparse representation of similarity/distance matrix 85 public double[] valP; // similarity/distance 86 public int[] rowP; // row index 87 public int[] colP; // col index 88 89 // for exact version: dense representation of distance/similarity matrix 90 public double[,] p; 91 92 // mapped data 93 public double[,] newData; 94 95 public int iter; 96 public double currentMomentum; 97 98 // helper variables (updated in each iteration) 99 public double[,] gains; 100 public double[,] uY; 101 public double[,] dY; 102 103 private TSNEState(TSNEState original, Cloner cloner) : base(original, cloner) { 104 } 105 public override IDeepCloneable Clone(Cloner cloner) { 106 return new TSNEState(this, cloner); 107 } 108 109 public TSNEState(T[] data, IDistance<T> distance, IRandom random, int newDimensions, double perplexity, double theta, int stopLyingIter, int momSwitchIter, double momentum, double finalMomentum, double eta) { 110 this.distance = distance; 111 this.random = random; 112 this.newDimensions = newDimensions; 113 this.perplexity = perplexity; 114 this.theta = theta; 115 this.stopLyingIter = stopLyingIter; 116 this.momSwitchIter = momSwitchIter; 117 this.currentMomentum = momentum; 118 this.finalMomentum = finalMomentum; 119 this.eta = eta; 120 121 122 // initialize 123 noDatapoints = data.Length; 124 if (noDatapoints  1 < 3 * perplexity) throw new ArgumentException("Perplexity too large for the number of data points!"); 125 126 exact = Math.Abs(theta) < double.Epsilon; 127 newData = new double[noDatapoints, newDimensions]; 128 dY = new double[noDatapoints, newDimensions]; 129 uY = new double[noDatapoints, newDimensions]; 130 gains = new double[noDatapoints, newDimensions]; 131 for (var i = 0; i < noDatapoints; i++) 132 for (var j = 0; j < newDimensions; j++) 133 gains[i, j] = 1.0; 134 135 p = null; 136 rowP = null; 137 colP = null; 138 valP = null; 139 140 //Calculate Similarities 141 if (exact) p = CalculateExactSimilarites(data, distance, perplexity); 142 else CalculateApproximateSimilarities(data, distance, perplexity, out rowP, out colP, out valP); 143 144 // Lie about the Pvalues 145 if (exact) for (var i = 0; i < noDatapoints; i++) for (var j = 0; j < noDatapoints; j++) p[i, j] *= 12.0; 146 else for (var i = 0; i < rowP[noDatapoints]; i++) valP[i] *= 12.0; 147 148 // Initialize solution (randomly) 149 var rand = new NormalDistributedRandom(random, 0, 1); 150 for (var i = 0; i < noDatapoints; i++) 151 for (var j = 0; j < newDimensions; j++) 152 newData[i, j] = rand.NextDouble() * .0001; // TODO const 153 } 154 155 public double EvaluateError() { 156 return exact ? EvaluateErrorExact(p, newData, noDatapoints, newDimensions) : EvaluateErrorApproximate(rowP, colP, valP, newData, theta); 157 } 158 159 private static void CalculateApproximateSimilarities(T[] data, IDistance<T> distance, double perplexity, out int[] rowP, out int[] colP, out double[] valP) { 160 // Compute asymmetric pairwise input similarities 161 ComputeGaussianPerplexity(data, distance, out rowP, out colP, out valP, perplexity, (int)(3 * perplexity)); // TODO: why 3? 162 // Symmetrize input similarities 163 int[] sRowP, symColP; 164 double[] sValP; 165 SymmetrizeMatrix(rowP, colP, valP, out sRowP, out symColP, out sValP); 166 rowP = sRowP; 167 colP = symColP; 168 valP = sValP; 169 var sumP = .0; 170 for (var i = 0; i < rowP[data.Length]; i++) sumP += valP[i]; 171 for (var i = 0; i < rowP[data.Length]; i++) valP[i] /= sumP; 172 } 173 private static double[,] CalculateExactSimilarites(T[] data, IDistance<T> distance, double perplexity) { 174 // Compute similarities 175 var p = new double[data.Length, data.Length]; 176 ComputeGaussianPerplexity(data, distance, p, perplexity); 177 // Symmetrize input similarities 178 for (var n = 0; n < data.Length; n++) { 179 for (var m = n + 1; m < data.Length; m++) { 180 p[n, m] += p[m, n]; 181 p[m, n] = p[n, m]; 182 } 183 } 184 var sumP = .0; 185 for (var i = 0; i < data.Length; i++) for (var j = 0; j < data.Length; j++) sumP += p[i, j]; 186 for (var i = 0; i < data.Length; i++) for (var j = 0; j < data.Length; j++) p[i, j] /= sumP; 187 return p; 188 } 189 190 private static void ComputeGaussianPerplexity(IReadOnlyList<T> x, IDistance<T> distance, out int[] rowP, out int[] colP, out double[] valP, double perplexity, int k) { 191 if (perplexity > k) throw new ArgumentException("Perplexity should be lower than k!"); 192 193 int n = x.Count; 194 // Allocate the memory we need 195 rowP = new int[n + 1]; 196 colP = new int[n * k]; 197 valP = new double[n * k]; 198 var curP = new double[n  1]; 199 rowP[0] = 0; 200 for (var i = 0; i < n; i++) rowP[i + 1] = rowP[i] + k; 201 202 var objX = new List<IndexedItem<T>>(); 203 for (var i = 0; i < n; i++) objX.Add(new IndexedItem<T>(i, x[i])); 204 205 // Build ball tree on data set 206 var tree = new VantagePointTree<IndexedItem<T>>(new IndexedItemDistance<T>(distance), objX); // do we really want to recreate the tree on each call? 207 208 // Loop over all points to find nearest neighbors 209 for (var i = 0; i < n; i++) { 210 IList<IndexedItem<T>> indices; 211 IList<double> distances; 212 213 // Find nearest neighbors 214 tree.Search(objX[i], k + 1, out indices, out distances); 215 216 // Initialize some variables for binary search 217 var found = false; 218 var beta = 1.0; 219 var minBeta = double.MinValue; 220 var maxBeta = double.MaxValue; 221 const double tol = 1e5; // TODO: why 1e5? 222 223 // Iterate until we found a good perplexity 224 var iter = 0; double sumP = 0; 225 while (!found && iter < 200) { 226 227 // Compute Gaussian kernel row 228 for (var m = 0; m < k; m++) curP[m] = Math.Exp(beta * distances[m + 1]); 229 230 // Compute entropy of current row 231 sumP = double.Epsilon; 232 for (var m = 0; m < k; m++) sumP += curP[m]; 233 var h = .0; 234 for (var m = 0; m < k; m++) h += beta * (distances[m + 1] * curP[m]); 235 h = h / sumP + Math.Log(sumP); 236 237 // Evaluate whether the entropy is within the tolerance level 238 var hdiff = h  Math.Log(perplexity); 239 if (hdiff < tol && hdiff < tol) { 240 found = true; 241 } else { 242 if (hdiff > 0) { 243 minBeta = beta; 244 if (maxBeta.IsAlmost(double.MaxValue)  maxBeta.IsAlmost(double.MinValue)) 245 beta *= 2.0; 246 else 247 beta = (beta + maxBeta) / 2.0; 248 } else { 249 maxBeta = beta; 250 if (minBeta.IsAlmost(double.MinValue)  minBeta.IsAlmost(double.MaxValue)) 251 beta /= 2.0; 252 else 253 beta = (beta + minBeta) / 2.0; 254 } 255 } 256 257 // Update iteration counter 258 iter++; 259 } 260 261 // Rownormalize current row of P and store in matrix 262 for (var m = 0; m < k; m++) curP[m] /= sumP; 263 for (var m = 0; m < k; m++) { 264 colP[rowP[i] + m] = indices[m + 1].Index; 265 valP[rowP[i] + m] = curP[m]; 266 } 267 } 268 } 269 private static void ComputeGaussianPerplexity(T[] x, IDistance<T> distance, double[,] p, double perplexity) { 270 // Compute the distance matrix 271 var dd = ComputeDistances(x, distance); 272 273 int n = x.Length; 274 // Compute the Gaussian kernel row by row 275 for (var i = 0; i < n; i++) { 276 // Initialize some variables 277 var found = false; 278 var beta = 1.0; 279 var minBeta = double.MaxValue; 280 var maxBeta = double.MaxValue; 281 const double tol = 1e5; 282 double sumP = 0; 283 284 // Iterate until we found a good perplexity 285 var iter = 0; 286 while (!found && iter < 200) { // TODO constant 287 288 // Compute Gaussian kernel row 289 for (var m = 0; m < n; m++) p[i, m] = Math.Exp(beta * dd[i][m]); 290 p[i, i] = double.Epsilon; 291 292 // Compute entropy of current row 293 sumP = double.Epsilon; 294 for (var m = 0; m < n; m++) sumP += p[i, m]; 295 var h = 0.0; 296 for (var m = 0; m < n; m++) h += beta * (dd[i][m] * p[i, m]); 297 h = h / sumP + Math.Log(sumP); 298 299 // Evaluate whether the entropy is within the tolerance level 300 var hdiff = h  Math.Log(perplexity); 301 if (hdiff < tol && hdiff < tol) { 302 found = true; 303 } else { 304 if (hdiff > 0) { 305 minBeta = beta; 306 if (maxBeta.IsAlmost(double.MaxValue)  maxBeta.IsAlmost(double.MinValue)) 307 beta *= 2.0; 308 else 309 beta = (beta + maxBeta) / 2.0; 310 } else { 311 maxBeta = beta; 312 if (minBeta.IsAlmost(double.MinValue)  minBeta.IsAlmost(double.MaxValue)) 313 beta /= 2.0; 314 else 315 beta = (beta + minBeta) / 2.0; 316 } 317 } 318 319 // Update iteration counter 320 iter++; 321 } 322 323 // Row normalize P 324 for (var m = 0; m < n; m++) p[i, m] /= sumP; 325 } 326 } 327 328 private static double[][] ComputeDistances(T[] x, IDistance<T> distance) { 329 return x.Select(m => x.Select(n => distance.Get(m, n)).ToArray()).ToArray(); 330 } 331 332 333 334 private static double EvaluateErrorExact(double[,] p, double[,] y, int n, int d) { 335 // Compute the squared Euclidean distance matrix 336 var dd = new double[n, n]; 337 var q = new double[n, n]; 338 ComputeSquaredEuclideanDistance(y, n, d, dd); 339 340 // Compute Qmatrix and normalization sum 341 var sumQ = double.Epsilon; 342 for (var n1 = 0; n1 < n; n1++) { 343 for (var m = 0; m < n; m++) { 344 if (n1 != m) { 345 q[n1, m] = 1 / (1 + dd[n1, m]); 346 sumQ += q[n1, m]; 347 } else q[n1, m] = double.Epsilon; 348 } 349 } 350 for (var i = 0; i < n; i++) for (var j = 0; j < n; j++) q[i, j] /= sumQ; 351 352 // Sum tSNE error 353 var c = .0; 354 for (var i = 0; i < n; i++) 355 for (var j = 0; j < n; j++) { 356 c += p[i, j] * Math.Log((p[i, j] + float.Epsilon) / (q[i, j] + float.Epsilon)); 357 } 358 return c; 359 } 360 private static double EvaluateErrorApproximate(IReadOnlyList<int> rowP, IReadOnlyList<int> colP, IReadOnlyList<double> valP, double[,] y, double theta) { 361 // Get estimate of normalization term 362 var n = y.GetLength(0); 363 var d = y.GetLength(1); 364 var tree = new SpacePartitioningTree(y); 365 var buff = new double[d]; 366 double sumQ = 0.0; 367 for (var i = 0; i < n; i++) tree.ComputeNonEdgeForces(i, theta, buff, ref sumQ); 368 369 // Loop over all edges to compute tSNE error 370 var c = .0; 371 for (var k = 0; k < n; k++) { 372 for (var i = rowP[k]; i < rowP[k + 1]; i++) { 373 var q = .0; 374 for (var j = 0; j < d; j++) buff[j] = y[k, j]; 375 for (var j = 0; j < d; j++) buff[j] = y[colP[i], j]; 376 for (var j = 0; j < d; j++) q += buff[j] * buff[j]; 377 q = 1.0 / (1.0 + q) / sumQ; 378 c += valP[i] * Math.Log((valP[i] + float.Epsilon) / (q + float.Epsilon)); 379 } 380 } 381 return c; 382 } 383 private static void SymmetrizeMatrix(IReadOnlyList<int> rowP, IReadOnlyList<int> colP, IReadOnlyList<double> valP, out int[] symRowP, out int[] symColP, out double[] symValP) { 384 385 // Count number of elements and row counts of symmetric matrix 386 var n = rowP.Count  1; 387 var rowCounts = new int[n]; 388 for (var j = 0; j < n; j++) { 389 for (var i = rowP[j]; i < rowP[j + 1]; i++) { 390 391 // Check whether element (col_P[i], n) is present 392 var present = false; 393 for (var m = rowP[colP[i]]; m < rowP[colP[i] + 1]; m++) { 394 if (colP[m] == j) present = true; 395 } 396 if (present) rowCounts[j]++; 397 else { 398 rowCounts[j]++; 399 rowCounts[colP[i]]++; 400 } 401 } 402 } 403 var noElem = 0; 404 for (var i = 0; i < n; i++) noElem += rowCounts[i]; 405 406 // Allocate memory for symmetrized matrix 407 symRowP = new int[n + 1]; 408 symColP = new int[noElem]; 409 symValP = new double[noElem]; 410 411 // Construct new row indices for symmetric matrix 412 symRowP[0] = 0; 413 for (var i = 0; i < n; i++) symRowP[i + 1] = symRowP[i] + rowCounts[i]; 414 415 // Fill the result matrix 416 var offset = new int[n]; 417 for (var j = 0; j < n; j++) { 418 for (var i = rowP[j]; i < rowP[j + 1]; i++) { // considering element(n, colP[i]) 419 420 // Check whether element (col_P[i], n) is present 421 var present = false; 422 for (var m = rowP[colP[i]]; m < rowP[colP[i] + 1]; m++) { 423 if (colP[m] != j) continue; 424 present = true; 425 if (j > colP[i]) continue; // make sure we do not add elements twice 426 symColP[symRowP[j] + offset[j]] = colP[i]; 427 symColP[symRowP[colP[i]] + offset[colP[i]]] = j; 428 symValP[symRowP[j] + offset[j]] = valP[i] + valP[m]; 429 symValP[symRowP[colP[i]] + offset[colP[i]]] = valP[i] + valP[m]; 430 } 431 432 // If (colP[i], n) is not present, there is no addition involved 433 if (!present) { 434 symColP[symRowP[j] + offset[j]] = colP[i]; 435 symColP[symRowP[colP[i]] + offset[colP[i]]] = j; 436 symValP[symRowP[j] + offset[j]] = valP[i]; 437 symValP[symRowP[colP[i]] + offset[colP[i]]] = valP[i]; 438 } 439 440 // Update offsets 441 if (present && (j > colP[i])) continue; 442 offset[j]++; 443 if (colP[i] != j) offset[colP[i]]++; 444 } 445 } 446 447 // Divide the result by two 448 for (var i = 0; i < noElem; i++) symValP[i] /= 2.0; 449 } 450 122 451 } 123 public override IDeepCloneable Clone(Cloner cloner) { return new TSNE<T>(this, cloner); } 124 public TSNE(IDistance<T> distance, IRandom random, ResultCollection results = null, int maxIter = 1000, int stopLyingIter = 250, int momSwitchIter = 250, double momentum = .5, double finalMomentum = .8, double eta = 200.0, Dictionary<string, List<int>> dataRowLookup = null, Dictionary<string, ScatterPlotDataRow> dataRows = null) { 125 this.distance = distance; 126 this.maxIter = maxIter; 127 this.stopLyingIter = stopLyingIter; 128 this.momSwitchIter = momSwitchIter; 129 this.momentum = momentum; 130 this.finalMomentum = finalMomentum; 131 this.eta = eta; 132 this.random = random; 133 this.results = results; 134 this.dataRowLookup = dataRowLookup; 135 if (dataRows != null) this.dataRows = dataRows; 136 else { this.dataRows = new Dictionary<string, ScatterPlotDataRow>(); } 452 453 public static TSNEState CreateState(T[] data, IDistance<T> distance, IRandom random, int newDimensions = 2, double perplexity = 25, double theta = 0, 454 int stopLyingIter = 250, int momSwitchIter = 250, double momentum = .5, double finalMomentum = .8, double eta = 200.0 455 ) { 456 return new TSNEState(data, distance, random, newDimensions, perplexity, theta, stopLyingIter, momSwitchIter, momentum, finalMomentum, eta); 137 457 } 138 #endregion 139 140 public double[,] Run(T[] data, int newDimensions, double perplexity, double theta) { 141 var currentMomentum = momentum; 142 var noDatapoints = data.Length; 143 if (noDatapoints  1 < 3 * perplexity) throw new ArgumentException("Perplexity too large for the number of data points!"); 144 SetUpResults(data); 145 Running = true; 146 var exact = Math.Abs(theta) < double.Epsilon; 147 var newData = new double[noDatapoints, newDimensions]; 148 var dY = new double[noDatapoints, newDimensions]; 149 var uY = new double[noDatapoints, newDimensions]; 150 var gains = new double[noDatapoints, newDimensions]; 151 for (var i = 0; i < noDatapoints; i++) for (var j = 0; j < newDimensions; j++) gains[i, j] = 1.0; 152 double[,] p = null; 153 int[] rowP = null; 154 int[] colP = null; 155 double[] valP = null; 156 var rand = new NormalDistributedRandom(random, 0, 1); 157 158 //Calculate Similarities 159 if (exact) p = CalculateExactSimilarites(data, perplexity); 160 else CalculateApproximateSimilarities(data, perplexity, out rowP, out colP, out valP); 161 162 // Lie about the Pvalues 163 if (exact) for (var i = 0; i < noDatapoints; i++) for (var j = 0; j < noDatapoints; j++) p[i, j] *= 12.0; 164 else for (var i = 0; i < rowP[noDatapoints]; i++) valP[i] *= 12.0; 165 166 // Initialize solution (randomly) 167 for (var i = 0; i < noDatapoints; i++) for (var j = 0; j < newDimensions; j++) newData[i, j] = rand.NextDouble() * .0001; // TODO const 168 169 // Perform main training loop 170 for (var iter = 0; iter < maxIter && Running; iter++) { 171 if (exact) ComputeExactGradient(p, newData, noDatapoints, newDimensions, dY); 172 else ComputeApproximateGradient(rowP, colP, valP, newData, noDatapoints, newDimensions, dY, theta); 173 // Update gains 174 for (var i = 0; i < noDatapoints; i++) for (var j = 0; j < newDimensions; j++) gains[i, j] = Math.Sign(dY[i, j]) != Math.Sign(uY[i, j]) ? gains[i, j] + .2 : gains[i, j] * .8; 175 for (var i = 0; i < noDatapoints; i++) for (var j = 0; j < newDimensions; j++) if (gains[i, j] < .01) gains[i, j] = .01; 176 // Perform gradient update (with momentum and gains) 177 for (var i = 0; i < noDatapoints; i++) for (var j = 0; j < newDimensions; j++) uY[i, j] = currentMomentum * uY[i, j]  eta * gains[i, j] * dY[i, j]; 178 for (var i = 0; i < noDatapoints; i++) for (var j = 0; j < newDimensions; j++) newData[i, j] = newData[i, j] + uY[i, j]; 179 // Make solution zeromean 180 ZeroMean(newData); 181 // Stop lying about the Pvalues after a while, and switch momentum 182 if (iter == stopLyingIter) { 183 if (exact) for (var i = 0; i < noDatapoints; i++) for (var j = 0; j < noDatapoints; j++) p[i, j] /= 12.0; 184 else for (var i = 0; i < rowP[noDatapoints]; i++) valP[i] /= 12.0; 185 } 186 if (iter == momSwitchIter) currentMomentum = finalMomentum; 187 188 Analyze(exact, iter, p, rowP, colP, valP, newData, noDatapoints, newDimensions, theta); 189 } 190 return newData; 458 459 460 public static double[,] Iterate(TSNEState state) { 461 if (state.exact) 462 ComputeExactGradient(state.p, state.newData, state.noDatapoints, state.newDimensions, state.dY); 463 else 464 ComputeApproximateGradient(state.rowP, state.colP, state.valP, state.newData, state.noDatapoints, state.newDimensions, state.dY, state.theta); 465 466 // Update gains 467 for (var i = 0; i < state.noDatapoints; i++) { 468 for (var j = 0; j < state.newDimensions; j++) { 469 state.gains[i, j] = Math.Sign(state.dY[i, j]) != Math.Sign(state.uY[i, j]) 470 ? state.gains[i, j] + .2 471 : state.gains[i, j] * .8; // 20% up or 20% down // TODO: +0.2?! 472 473 if (state.gains[i, j] < .01) state.gains[i, j] = .01; 474 } 475 } 476 477 478 // Perform gradient update (with momentum and gains) 479 for (var i = 0; i < state.noDatapoints; i++) 480 for (var j = 0; j < state.newDimensions; j++) 481 state.uY[i, j] = state.currentMomentum * state.uY[i, j]  state.eta * state.gains[i, j] * state.dY[i, j]; 482 483 for (var i = 0; i < state.noDatapoints; i++) 484 for (var j = 0; j < state.newDimensions; j++) 485 state.newData[i, j] = state.newData[i, j] + state.uY[i, j]; 486 487 // Make solution zeromean 488 ZeroMean(state.newData); 489 // Stop lying about the Pvalues after a while, and switch momentum 490 491 if (state.iter == state.stopLyingIter) { 492 if (state.exact) 493 for (var i = 0; i < state.noDatapoints; i++) for (var j = 0; j < state.noDatapoints; j++) state.p[i, j] /= 12.0; //XXX why 12? 494 else 495 for (var i = 0; i < state.rowP[state.noDatapoints]; i++) state.valP[i] /= 12.0; // XXX are we not scaling all values? 496 } 497 498 if (state.iter == state.momSwitchIter) 499 state.currentMomentum = state.finalMomentum; 500 501 state.iter++; 502 return state.newData; 191 503 } 192 504 193 #region helpers 194 195 private void SetUpResults(IReadOnlyCollection<T> data) { 196 if (dataRowLookup == null) dataRowLookup = new Dictionary<string, List<int>> { { "Data", Enumerable.Range(0, data.Count).ToList() } }; 197 if (results == null) return; 198 199 if (!results.ContainsKey(IterationResultName)) results.Add(new Result(IterationResultName, new IntValue(0))); 200 else ((IntValue)results[IterationResultName].Value).Value = 0; 201 202 if (!results.ContainsKey(ErrorResultName)) results.Add(new Result(ErrorResultName, new DoubleValue(0))); 203 else ((DoubleValue)results[ErrorResultName].Value).Value = 0; 204 205 if (!results.ContainsKey(ErrorPlotResultName)) results.Add(new Result(ErrorPlotResultName, new DataTable(ErrorPlotResultName, "Development of errors during gradient descent"))); 206 else results[ErrorPlotResultName].Value = new DataTable(ErrorPlotResultName, "Development of errors during gradient descent"); 207 208 var plot = results[ErrorPlotResultName].Value as DataTable; 209 if (plot == null) throw new ArgumentException("could not create/access error data table in results collection"); 210 211 if (!plot.Rows.ContainsKey("errors")) plot.Rows.Add(new DataRow("errors")); 212 plot.Rows["errors"].Values.Clear(); 213 214 results.Add(new Result(ScatterPlotResultName, "Plot of the projected data", new ScatterPlot(DataResultName, ""))); 215 results.Add(new Result(DataResultName, "Projected Data", new DoubleMatrix())); 216 217 } 218 private void Analyze(bool exact, int iter, double[,] p, int[] rowP, int[] colP, double[] valP, double[,] newData, int noDatapoints, int newDimensions, double theta) { 219 if (results == null) return; 220 var plot = results[ErrorPlotResultName].Value as DataTable; 221 if (plot == null) throw new ArgumentException("Could not create/access error data table in results collection. Was it removed by some effect?"); 222 var errors = plot.Rows["errors"].Values; 223 var c = exact 224 ? EvaluateErrorExact(p, newData, noDatapoints, newDimensions) 225 : EvaluateErrorApproximate(rowP, colP, valP, newData, theta); 226 errors.Add(c); 227 ((IntValue)results[IterationResultName].Value).Value = iter + 1; 228 ((DoubleValue)results[ErrorResultName].Value).Value = errors.Last(); 229 230 var ndata = Normalize(newData); 231 results[DataResultName].Value = new DoubleMatrix(ndata); 232 var splot = results[ScatterPlotResultName].Value as ScatterPlot; 233 FillScatterPlot(ndata, splot); 234 235 236 } 237 private void FillScatterPlot(double[,] lowDimData, ScatterPlot plot) { 238 foreach (var rowName in dataRowLookup.Keys) { 239 if (!plot.Rows.ContainsKey(rowName)) 240 plot.Rows.Add(dataRows.ContainsKey(rowName) ? dataRows[rowName] : new ScatterPlotDataRow(rowName, "", new List<Point2D<double>>())); 241 plot.Rows[rowName].Points.Replace(dataRowLookup[rowName].Select(i => new Point2D<double>(lowDimData[i, 0], lowDimData[i, 1]))); 242 } 243 } 244 private static double[,] Normalize(double[,] data) { 245 var max = new double[data.GetLength(1)]; 246 var min = new double[data.GetLength(1)]; 247 var res = new double[data.GetLength(0), data.GetLength(1)]; 248 for (var i = 0; i < max.Length; i++) max[i] = min[i] = data[0, i]; 249 for (var i = 0; i < data.GetLength(0); i++) 250 for (var j = 0; j < data.GetLength(1); j++) { 251 var v = data[i, j]; 252 max[j] = Math.Max(max[j], v); 253 min[j] = Math.Min(min[j], v); 254 } 255 for (var i = 0; i < data.GetLength(0); i++) { 256 for (var j = 0; j < data.GetLength(1); j++) { 257 res[i, j] = (data[i, j]  (max[j] + min[j]) / 2) / (max[j]  min[j]); 258 } 259 } 260 return res; 261 } 262 private void CalculateApproximateSimilarities(T[] data, double perplexity, out int[] rowP, out int[] colP, out double[] valP) { 263 // Compute asymmetric pairwise input similarities 264 ComputeGaussianPerplexity(data, data.Length, out rowP, out colP, out valP, perplexity, (int)(3 * perplexity)); 265 // Symmetrize input similarities 266 int[] sRowP, symColP; 267 double[] sValP; 268 SymmetrizeMatrix(rowP, colP, valP, out sRowP, out symColP, out sValP); 269 rowP = sRowP; 270 colP = symColP; 271 valP = sValP; 272 var sumP = .0; 273 for (var i = 0; i < rowP[data.Length]; i++) sumP += valP[i]; 274 for (var i = 0; i < rowP[data.Length]; i++) valP[i] /= sumP; 275 } 276 private double[,] CalculateExactSimilarites(T[] data, double perplexity) { 277 // Compute similarities 278 var p = new double[data.Length, data.Length]; 279 ComputeGaussianPerplexity(data, data.Length, p, perplexity); 280 // Symmetrize input similarities 281 for (var n = 0; n < data.Length; n++) { 282 for (var m = n + 1; m < data.Length; m++) { 283 p[n, m] += p[m, n]; 284 p[m, n] = p[n, m]; 285 } 286 } 287 var sumP = .0; 288 for (var i = 0; i < data.Length; i++) for (var j = 0; j < data.Length; j++) sumP += p[i, j]; 289 for (var i = 0; i < data.Length; i++) for (var j = 0; j < data.Length; j++) p[i, j] /= sumP; 290 return p; 291 } 292 293 private void ComputeGaussianPerplexity(IReadOnlyList<T> x, int n, out int[] rowP, out int[] colP, out double[] valP, double perplexity, int k) { 294 if (perplexity > k) throw new ArgumentException("Perplexity should be lower than K!"); 295 296 // Allocate the memory we need 297 rowP = new int[n + 1]; 298 colP = new int[n * k]; 299 valP = new double[n * k]; 300 var curP = new double[n  1]; 301 rowP[0] = 0; 302 for (var i = 0; i < n; i++) rowP[i + 1] = rowP[i] + k; 303 304 var objX = new List<IndexedItem<T>>(); 305 for (var i = 0; i < n; i++) objX.Add(new IndexedItem<T>(i, x[i])); 306 307 // Build ball tree on data set 308 var tree = new VantagePointTree<IndexedItem<T>>(new IndexedItemDistance<T>(distance), objX); // do we really want to recreate the tree on each call? 309 310 // Loop over all points to find nearest neighbors 311 for (var i = 0; i < n; i++) { 312 IList<IndexedItem<T>> indices; 313 IList<double> distances; 314 315 // Find nearest neighbors 316 tree.Search(objX[i], k + 1, out indices, out distances); 317 318 // Initialize some variables for binary search 319 var found = false; 320 var beta = 1.0; 321 var minBeta = double.MinValue; 322 var maxBeta = double.MaxValue; 323 const double tol = 1e5; // TODO: why 1e5? 324 325 // Iterate until we found a good perplexity 326 var iter = 0; double sumP = 0; 327 while (!found && iter < 200) { 328 329 // Compute Gaussian kernel row 330 for (var m = 0; m < k; m++) curP[m] = Math.Exp(beta * distances[m + 1]); 331 332 // Compute entropy of current row 333 sumP = double.Epsilon; 334 for (var m = 0; m < k; m++) sumP += curP[m]; 335 var h = .0; 336 for (var m = 0; m < k; m++) h += beta * (distances[m + 1] * curP[m]); 337 h = h / sumP + Math.Log(sumP); 338 339 // Evaluate whether the entropy is within the tolerance level 340 var hdiff = h  Math.Log(perplexity); 341 if (hdiff < tol && hdiff < tol) { 342 found = true; 343 } else { 344 if (hdiff > 0) { 345 minBeta = beta; 346 if (maxBeta.IsAlmost(double.MaxValue)  maxBeta.IsAlmost(double.MinValue)) 347 beta *= 2.0; 348 else 349 beta = (beta + maxBeta) / 2.0; 350 } else { 351 maxBeta = beta; 352 if (minBeta.IsAlmost(double.MinValue)  minBeta.IsAlmost(double.MaxValue)) 353 beta /= 2.0; 354 else 355 beta = (beta + minBeta) / 2.0; 356 } 357 } 358 359 // Update iteration counter 360 iter++; 361 } 362 363 // Rownormalize current row of P and store in matrix 364 for (var m = 0; m < k; m++) curP[m] /= sumP; 365 for (var m = 0; m < k; m++) { 366 colP[rowP[i] + m] = indices[m + 1].Index; 367 valP[rowP[i] + m] = curP[m]; 368 } 369 } 370 } 371 private void ComputeGaussianPerplexity(T[] x, int n, double[,] p, double perplexity) { 372 // Compute the distance matrix 373 var dd = ComputeDistances(x); 374 375 // Compute the Gaussian kernel row by row 376 for (var i = 0; i < n; i++) { 377 // Initialize some variables 378 var found = false; 379 var beta = 1.0; 380 var minBeta = double.MaxValue; 381 var maxBeta = double.MaxValue; 382 const double tol = 1e5; 383 double sumP = 0; 384 385 // Iterate until we found a good perplexity 386 var iter = 0; 387 while (!found && iter < 200) { // TODO constant 388 389 // Compute Gaussian kernel row 390 for (var m = 0; m < n; m++) p[i, m] = Math.Exp(beta * dd[i][m]); 391 p[i, i] = double.Epsilon; 392 393 // Compute entropy of current row 394 sumP = double.Epsilon; 395 for (var m = 0; m < n; m++) sumP += p[i, m]; 396 var h = 0.0; 397 for (var m = 0; m < n; m++) h += beta * (dd[i][m] * p[i, m]); 398 h = h / sumP + Math.Log(sumP); 399 400 // Evaluate whether the entropy is within the tolerance level 401 var hdiff = h  Math.Log(perplexity); 402 if (hdiff < tol && hdiff < tol) { 403 found = true; 404 } else { 405 if (hdiff > 0) { 406 minBeta = beta; 407 if (maxBeta.IsAlmost(double.MaxValue)  maxBeta.IsAlmost(double.MinValue)) 408 beta *= 2.0; 409 else 410 beta = (beta + maxBeta) / 2.0; 411 } else { 412 maxBeta = beta; 413 if (minBeta.IsAlmost(double.MinValue)  minBeta.IsAlmost(double.MaxValue)) 414 beta /= 2.0; 415 else 416 beta = (beta + minBeta) / 2.0; 417 } 418 } 419 420 // Update iteration counter 421 iter++; 422 } 423 424 // Row normalize P 425 for (var m = 0; m < n; m++) p[i, m] /= sumP; 426 } 427 } 428 429 private double[][] ComputeDistances(T[] x) { 430 return x.Select(m => x.Select(n => distance.Get(m, n)).ToArray()).ToArray(); 505 506 private static void ComputeApproximateGradient(int[] rowP, int[] colP, double[] valP, double[,] y, int n, int d, double[,] dC, double theta) { 507 var tree = new SpacePartitioningTree(y); 508 double sumQ = 0.0; 509 var posF = new double[n, d]; 510 var negF = new double[n, d]; 511 tree.ComputeEdgeForces(rowP, colP, valP, n, posF); 512 var row = new double[d]; 513 for (var n1 = 0; n1 < n; n1++) { 514 Buffer.BlockCopy(negF, (sizeof(double) * n1 * d), row, 0, d); 515 tree.ComputeNonEdgeForces(n1, theta, row, ref sumQ); 516 } 517 518 // Compute final tSNE gradient 519 for (var i = 0; i < n; i++) 520 for (var j = 0; j < d; j++) { 521 dC[i, j] = posF[i, j]  negF[i, j] / sumQ; 522 } 431 523 } 432 524 … … 462 554 } 463 555 } 556 464 557 private static void ComputeSquaredEuclideanDistance(double[,] x, int n, int d, double[,] dd) { 465 558 var dataSums = new double[n]; … … 485 578 } 486 579 } 487 private static void ComputeApproximateGradient(int[] rowP, int[] colP, double[] valP, double[,] y, int n, int d, double[,] dC, double theta) { 488 var tree = new SpacePartitioningTree(y); 489 double[] sumQ = { 0 }; 490 var posF = new double[n, d]; 491 var negF = new double[n, d]; 492 tree.ComputeEdgeForces(rowP, colP, valP, n, posF); 493 var row = new double[d]; 494 for (var n1 = 0; n1 < n; n1++) { 495 Buffer.BlockCopy(negF, (sizeof(double) * n1 * d), row, 0, d); 496 tree.ComputeNonEdgeForces(n1, theta, row, sumQ); 497 } 498 499 // Compute final tSNE gradient 500 for (var i = 0; i < n; i++) 501 for (var j = 0; j < d; j++) { 502 dC[i, j] = (posF[i, j]  negF[i, j]) / sumQ[0]; // TODO: check parenthesis 503 } 504 } 505 506 private static double EvaluateErrorExact(double[,] p, double[,] y, int n, int d) { 507 // Compute the squared Euclidean distance matrix 508 var dd = new double[n, n]; 509 var q = new double[n, n]; 510 ComputeSquaredEuclideanDistance(y, n, d, dd); 511 512 // Compute Qmatrix and normalization sum 513 var sumQ = double.Epsilon; 514 for (var n1 = 0; n1 < n; n1++) { 515 for (var m = 0; m < n; m++) { 516 if (n1 != m) { 517 q[n1, m] = 1 / (1 + dd[n1, m]); 518 sumQ += q[n1, m]; 519 } else q[n1, m] = double.Epsilon; 520 } 521 } 522 for (var i = 0; i < n; i++) for (var j = 0; j < n; j++) q[i, j] /= sumQ; 523 524 // Sum tSNE error 525 var c = .0; 526 for (var i = 0; i < n; i++) 527 for (var j = 0; j < n; j++) { 528 c += p[i, j] * Math.Log((p[i, j] + float.Epsilon) / (q[i, j] + float.Epsilon)); 529 } 530 return c; 531 } 532 private static double EvaluateErrorApproximate(IReadOnlyList<int> rowP, IReadOnlyList<int> colP, IReadOnlyList<double> valP, double[,] y, double theta) { 533 // Get estimate of normalization term 534 var n = y.GetLength(0); 535 var d = y.GetLength(1); 536 var tree = new SpacePartitioningTree(y); 537 var buff = new double[d]; 538 double[] sumQ = { 0 }; 539 for (var i = 0; i < n; i++) tree.ComputeNonEdgeForces(i, theta, buff, sumQ); 540 541 // Loop over all edges to compute tSNE error 542 var c = .0; 543 for (var k = 0; k < n; k++) { 544 for (var i = rowP[k]; i < rowP[k + 1]; i++) { 545 var q = .0; 546 for (var j = 0; j < d; j++) buff[j] = y[k, j]; 547 for (var j = 0; j < d; j++) buff[j] = y[colP[i], j]; 548 for (var j = 0; j < d; j++) q += buff[j] * buff[j]; 549 q = 1.0 / (1.0 + q) / sumQ[0]; 550 c += valP[i] * Math.Log((valP[i] + float.Epsilon) / (q + float.Epsilon)); 551 } 552 } 553 return c; 554 } 555 private static void SymmetrizeMatrix(IReadOnlyList<int> rowP, IReadOnlyList<int> colP, IReadOnlyList<double> valP, out int[] symRowP, out int[] symColP, out double[] symValP) { 556 557 // Count number of elements and row counts of symmetric matrix 558 var n = rowP.Count  1; 559 var rowCounts = new int[n]; 560 for (var j = 0; j < n; j++) { 561 for (var i = rowP[j]; i < rowP[j + 1]; i++) { 562 563 // Check whether element (col_P[i], n) is present 564 var present = false; 565 for (var m = rowP[colP[i]]; m < rowP[colP[i] + 1]; m++) { 566 if (colP[m] == j) present = true; 567 } 568 if (present) rowCounts[j]++; 569 else { 570 rowCounts[j]++; 571 rowCounts[colP[i]]++; 572 } 573 } 574 } 575 var noElem = 0; 576 for (var i = 0; i < n; i++) noElem += rowCounts[i]; 577 578 // Allocate memory for symmetrized matrix 579 symRowP = new int[n + 1]; 580 symColP = new int[noElem]; 581 symValP = new double[noElem]; 582 583 // Construct new row indices for symmetric matrix 584 symRowP[0] = 0; 585 for (var i = 0; i < n; i++) symRowP[i + 1] = symRowP[i] + rowCounts[i]; 586 587 // Fill the result matrix 588 var offset = new int[n]; 589 for (var j = 0; j < n; j++) { 590 for (var i = rowP[j]; i < rowP[j + 1]; i++) { // considering element(n, colP[i]) 591 592 // Check whether element (col_P[i], n) is present 593 var present = false; 594 for (var m = rowP[colP[i]]; m < rowP[colP[i] + 1]; m++) { 595 if (colP[m] != j) continue; 596 present = true; 597 if (j > colP[i]) continue; // make sure we do not add elements twice 598 symColP[symRowP[j] + offset[j]] = colP[i]; 599 symColP[symRowP[colP[i]] + offset[colP[i]]] = j; 600 symValP[symRowP[j] + offset[j]] = valP[i] + valP[m]; 601 symValP[symRowP[colP[i]] + offset[colP[i]]] = valP[i] + valP[m]; 602 } 603 604 // If (colP[i], n) is not present, there is no addition involved 605 if (!present) { 606 symColP[symRowP[j] + offset[j]] = colP[i]; 607 symColP[symRowP[colP[i]] + offset[colP[i]]] = j; 608 symValP[symRowP[j] + offset[j]] = valP[i]; 609 symValP[symRowP[colP[i]] + offset[colP[i]]] = valP[i]; 610 } 611 612 // Update offsets 613 if (present && (j > colP[i])) continue; 614 offset[j]++; 615 if (colP[i] != j) offset[colP[i]]++; 616 } 617 } 618 619 // Divide the result by two 620 for (var i = 0; i < noElem; i++) symValP[i] /= 2.0; 621 } 580 622 581 private static void ZeroMean(double[,] x) { 623 582 // Compute data mean … … 640 599 } 641 600 } 642 #endregion643 601 } 644 602 }
Note: See TracChangeset
for help on using the changeset viewer.