[14414] | 1 | #region License Information
|
---|
| 2 | /* HeuristicLab
|
---|
| 3 | * Copyright (C) 2002-2016 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 |
|
---|
| 21 | //Code is based on an implementation from Laurens van der Maaten
|
---|
| 22 |
|
---|
| 23 | /*
|
---|
| 24 | *
|
---|
| 25 | * Copyright (c) 2014, Laurens van der Maaten (Delft University of Technology)
|
---|
| 26 | * All rights reserved.
|
---|
| 27 | *
|
---|
| 28 | * Redistribution and use in source and binary forms, with or without
|
---|
| 29 | * modification, are permitted provided that the following conditions are met:
|
---|
| 30 | * 1. Redistributions of source code must retain the above copyright
|
---|
| 31 | * notice, this list of conditions and the following disclaimer.
|
---|
| 32 | * 2. Redistributions in binary form must reproduce the above copyright
|
---|
| 33 | * notice, this list of conditions and the following disclaimer in the
|
---|
| 34 | * documentation and/or other materials provided with the distribution.
|
---|
| 35 | * 3. All advertising materials mentioning features or use of this software
|
---|
| 36 | * must display the following acknowledgement:
|
---|
| 37 | * This product includes software developed by the Delft University of Technology.
|
---|
| 38 | * 4. Neither the name of the Delft University of Technology nor the names of
|
---|
| 39 | * its contributors may be used to endorse or promote products derived from
|
---|
| 40 | * this software without specific prior written permission.
|
---|
| 41 | *
|
---|
| 42 | * THIS SOFTWARE IS PROVIDED BY LAURENS VAN DER MAATEN ''AS IS'' AND ANY EXPRESS
|
---|
| 43 | * OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
|
---|
| 44 | * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO
|
---|
| 45 | * EVENT SHALL LAURENS VAN DER MAATEN BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
|
---|
| 46 | * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
|
---|
| 47 | * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR
|
---|
| 48 | * BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
|
---|
| 49 | * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING
|
---|
| 50 | * IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY
|
---|
| 51 | * OF SUCH DAMAGE.
|
---|
| 52 | *
|
---|
| 53 | */
|
---|
| 54 | #endregion
|
---|
| 55 |
|
---|
| 56 | using System;
|
---|
| 57 | using System.Collections.Generic;
|
---|
| 58 | using System.Linq;
|
---|
[14785] | 59 | using HeuristicLab.Collections;
|
---|
[14414] | 60 | using HeuristicLab.Common;
|
---|
| 61 | using HeuristicLab.Core;
|
---|
| 62 | using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
|
---|
| 63 | using HeuristicLab.Random;
|
---|
| 64 |
|
---|
| 65 | namespace HeuristicLab.Algorithms.DataAnalysis {
|
---|
| 66 | [StorableClass]
|
---|
[14788] | 67 | public class TSNE<T> {
|
---|
[14414] | 68 |
|
---|
[14788] | 69 | [StorableClass]
|
---|
| 70 | public sealed class TSNEState : DeepCloneable {
|
---|
| 71 | // initialized once
|
---|
[14806] | 72 | [Storable]
|
---|
[14788] | 73 | public IDistance<T> distance;
|
---|
[14806] | 74 | [Storable]
|
---|
[14788] | 75 | public IRandom random;
|
---|
[14806] | 76 | [Storable]
|
---|
[14788] | 77 | public double perplexity;
|
---|
[14806] | 78 | [Storable]
|
---|
[14788] | 79 | public bool exact;
|
---|
[14806] | 80 | [Storable]
|
---|
[14788] | 81 | public int noDatapoints;
|
---|
[14806] | 82 | [Storable]
|
---|
[14788] | 83 | public double finalMomentum;
|
---|
[14806] | 84 | [Storable]
|
---|
[14788] | 85 | public int momSwitchIter;
|
---|
[14806] | 86 | [Storable]
|
---|
[14788] | 87 | public int stopLyingIter;
|
---|
[14806] | 88 | [Storable]
|
---|
[14788] | 89 | public double theta;
|
---|
[14806] | 90 | [Storable]
|
---|
[14788] | 91 | public double eta;
|
---|
[14806] | 92 | [Storable]
|
---|
[14788] | 93 | public int newDimensions;
|
---|
[14414] | 94 |
|
---|
[14788] | 95 | // for approximate version: sparse representation of similarity/distance matrix
|
---|
[14806] | 96 | [Storable]
|
---|
[14788] | 97 | public double[] valP; // similarity/distance
|
---|
[14806] | 98 | [Storable]
|
---|
[14788] | 99 | public int[] rowP; // row index
|
---|
[14806] | 100 | [Storable]
|
---|
[14788] | 101 | public int[] colP; // col index
|
---|
[14414] | 102 |
|
---|
[14788] | 103 | // for exact version: dense representation of distance/similarity matrix
|
---|
[14806] | 104 | [Storable]
|
---|
[14788] | 105 | public double[,] p;
|
---|
[14512] | 106 |
|
---|
[14788] | 107 | // mapped data
|
---|
[14806] | 108 | [Storable]
|
---|
[14788] | 109 | public double[,] newData;
|
---|
[14414] | 110 |
|
---|
[14806] | 111 | [Storable]
|
---|
[14788] | 112 | public int iter;
|
---|
[14806] | 113 | [Storable]
|
---|
[14788] | 114 | public double currentMomentum;
|
---|
[14414] | 115 |
|
---|
[14788] | 116 | // helper variables (updated in each iteration)
|
---|
[14806] | 117 | [Storable]
|
---|
[14788] | 118 | public double[,] gains;
|
---|
[14806] | 119 | [Storable]
|
---|
[14788] | 120 | public double[,] uY;
|
---|
[14806] | 121 | [Storable]
|
---|
[14788] | 122 | public double[,] dY;
|
---|
[14512] | 123 |
|
---|
[14788] | 124 | private TSNEState(TSNEState original, Cloner cloner) : base(original, cloner) {
|
---|
[14806] | 125 | this.distance = cloner.Clone(original.distance);
|
---|
| 126 | this.random = cloner.Clone(original.random);
|
---|
| 127 | this.perplexity = original.perplexity;
|
---|
| 128 | this.exact = original.exact;
|
---|
| 129 | this.noDatapoints = original.noDatapoints;
|
---|
| 130 | this.finalMomentum = original.finalMomentum;
|
---|
| 131 | this.momSwitchIter = original.momSwitchIter;
|
---|
| 132 | this.stopLyingIter = original.stopLyingIter;
|
---|
| 133 | this.theta = original.theta;
|
---|
| 134 | this.eta = original.eta;
|
---|
| 135 | this.newDimensions = original.newDimensions;
|
---|
| 136 | if(original.valP != null) {
|
---|
| 137 | this.valP = new double[original.valP.Length];
|
---|
| 138 | Array.Copy(original.valP, this.valP, this.valP.Length);
|
---|
| 139 | }
|
---|
| 140 | if(original.rowP != null) {
|
---|
| 141 | this.rowP = new int[original.rowP.Length];
|
---|
| 142 | Array.Copy(original.rowP, this.rowP, this.rowP.Length);
|
---|
| 143 | }
|
---|
| 144 | if(original.colP != null) {
|
---|
| 145 | this.colP = new int[original.colP.Length];
|
---|
| 146 | Array.Copy(original.colP, this.colP, this.colP.Length);
|
---|
| 147 | }
|
---|
| 148 | if(original.p != null) {
|
---|
| 149 | this.p = new double[original.p.GetLength(0), original.p.GetLength(1)];
|
---|
| 150 | Array.Copy(original.p, this.p, this.p.Length);
|
---|
| 151 | }
|
---|
| 152 | this.newData = new double[original.newData.GetLength(0), original.newData.GetLength(1)];
|
---|
| 153 | Array.Copy(original.newData, this.newData, this.newData.Length);
|
---|
| 154 | this.iter = original.iter;
|
---|
| 155 | this.currentMomentum = original.currentMomentum;
|
---|
| 156 | this.gains = new double[original.gains.GetLength(0), original.gains.GetLength(1)];
|
---|
| 157 | Array.Copy(original.gains, this.gains, this.gains.Length);
|
---|
| 158 | this.uY = new double[original.uY.GetLength(0), original.uY.GetLength(1)];
|
---|
| 159 | Array.Copy(original.uY, this.uY, this.uY.Length);
|
---|
| 160 | this.dY = new double[original.dY.GetLength(0), original.dY.GetLength(1)];
|
---|
| 161 | Array.Copy(original.dY, this.dY, this.dY.Length);
|
---|
[14788] | 162 | }
|
---|
[14806] | 163 |
|
---|
[14788] | 164 | public override IDeepCloneable Clone(Cloner cloner) {
|
---|
| 165 | return new TSNEState(this, cloner);
|
---|
| 166 | }
|
---|
[14414] | 167 |
|
---|
[14788] | 168 | 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) {
|
---|
| 169 | this.distance = distance;
|
---|
| 170 | this.random = random;
|
---|
| 171 | this.newDimensions = newDimensions;
|
---|
| 172 | this.perplexity = perplexity;
|
---|
| 173 | this.theta = theta;
|
---|
| 174 | this.stopLyingIter = stopLyingIter;
|
---|
| 175 | this.momSwitchIter = momSwitchIter;
|
---|
| 176 | this.currentMomentum = momentum;
|
---|
| 177 | this.finalMomentum = finalMomentum;
|
---|
| 178 | this.eta = eta;
|
---|
[14414] | 179 |
|
---|
| 180 |
|
---|
[14788] | 181 | // initialize
|
---|
| 182 | noDatapoints = data.Length;
|
---|
[14806] | 183 | if(noDatapoints - 1 < 3 * perplexity)
|
---|
| 184 | throw new ArgumentException("Perplexity too large for the number of data points!");
|
---|
[14788] | 185 |
|
---|
| 186 | exact = Math.Abs(theta) < double.Epsilon;
|
---|
| 187 | newData = new double[noDatapoints, newDimensions];
|
---|
| 188 | dY = new double[noDatapoints, newDimensions];
|
---|
| 189 | uY = new double[noDatapoints, newDimensions];
|
---|
| 190 | gains = new double[noDatapoints, newDimensions];
|
---|
[14806] | 191 | for(var i = 0; i < noDatapoints; i++)
|
---|
| 192 | for(var j = 0; j < newDimensions; j++)
|
---|
[14788] | 193 | gains[i, j] = 1.0;
|
---|
| 194 |
|
---|
| 195 | p = null;
|
---|
| 196 | rowP = null;
|
---|
| 197 | colP = null;
|
---|
| 198 | valP = null;
|
---|
| 199 |
|
---|
| 200 | //Calculate Similarities
|
---|
[14806] | 201 | if(exact) p = CalculateExactSimilarites(data, distance, perplexity);
|
---|
[14788] | 202 | else CalculateApproximateSimilarities(data, distance, perplexity, out rowP, out colP, out valP);
|
---|
| 203 |
|
---|
| 204 | // Lie about the P-values
|
---|
[14806] | 205 | if(exact) for(var i = 0; i < noDatapoints; i++) for(var j = 0; j < noDatapoints; j++) p[i, j] *= 12.0;
|
---|
| 206 | else for(var i = 0; i < rowP[noDatapoints]; i++) valP[i] *= 12.0;
|
---|
[14788] | 207 |
|
---|
| 208 | // Initialize solution (randomly)
|
---|
| 209 | var rand = new NormalDistributedRandom(random, 0, 1);
|
---|
[14806] | 210 | for(var i = 0; i < noDatapoints; i++)
|
---|
| 211 | for(var j = 0; j < newDimensions; j++)
|
---|
| 212 | newData[i, j] = rand.NextDouble() * .0001; // TODO const ?
|
---|
[14414] | 213 | }
|
---|
| 214 |
|
---|
[14788] | 215 | public double EvaluateError() {
|
---|
[14806] | 216 | return exact ?
|
---|
| 217 | EvaluateErrorExact(p, newData, noDatapoints, newDimensions) :
|
---|
| 218 | EvaluateErrorApproximate(rowP, colP, valP, newData, theta);
|
---|
[14788] | 219 | }
|
---|
[14512] | 220 |
|
---|
[14788] | 221 | private static void CalculateApproximateSimilarities(T[] data, IDistance<T> distance, double perplexity, out int[] rowP, out int[] colP, out double[] valP) {
|
---|
| 222 | // Compute asymmetric pairwise input similarities
|
---|
| 223 | ComputeGaussianPerplexity(data, distance, out rowP, out colP, out valP, perplexity, (int)(3 * perplexity)); // TODO: why 3?
|
---|
| 224 | // Symmetrize input similarities
|
---|
| 225 | int[] sRowP, symColP;
|
---|
| 226 | double[] sValP;
|
---|
| 227 | SymmetrizeMatrix(rowP, colP, valP, out sRowP, out symColP, out sValP);
|
---|
| 228 | rowP = sRowP;
|
---|
| 229 | colP = symColP;
|
---|
| 230 | valP = sValP;
|
---|
| 231 | var sumP = .0;
|
---|
[14806] | 232 | for(var i = 0; i < rowP[data.Length]; i++) sumP += valP[i];
|
---|
| 233 | for(var i = 0; i < rowP[data.Length]; i++) valP[i] /= sumP;
|
---|
[14788] | 234 | }
|
---|
[14806] | 235 |
|
---|
[14788] | 236 | private static double[,] CalculateExactSimilarites(T[] data, IDistance<T> distance, double perplexity) {
|
---|
| 237 | // Compute similarities
|
---|
| 238 | var p = new double[data.Length, data.Length];
|
---|
| 239 | ComputeGaussianPerplexity(data, distance, p, perplexity);
|
---|
| 240 | // Symmetrize input similarities
|
---|
[14806] | 241 | for(var n = 0; n < data.Length; n++) {
|
---|
| 242 | for(var m = n + 1; m < data.Length; m++) {
|
---|
[14788] | 243 | p[n, m] += p[m, n];
|
---|
| 244 | p[m, n] = p[n, m];
|
---|
| 245 | }
|
---|
| 246 | }
|
---|
| 247 | var sumP = .0;
|
---|
[14806] | 248 | for(var i = 0; i < data.Length; i++) for(var j = 0; j < data.Length; j++) sumP += p[i, j];
|
---|
| 249 | for(var i = 0; i < data.Length; i++) for(var j = 0; j < data.Length; j++) p[i, j] /= sumP;
|
---|
[14788] | 250 | return p;
|
---|
| 251 | }
|
---|
[14742] | 252 |
|
---|
[14788] | 253 | private static void ComputeGaussianPerplexity(IReadOnlyList<T> x, IDistance<T> distance, out int[] rowP, out int[] colP, out double[] valP, double perplexity, int k) {
|
---|
[14806] | 254 | if(perplexity > k) throw new ArgumentException("Perplexity should be lower than k!");
|
---|
[14512] | 255 |
|
---|
[14788] | 256 | int n = x.Count;
|
---|
| 257 | // Allocate the memory we need
|
---|
| 258 | rowP = new int[n + 1];
|
---|
| 259 | colP = new int[n * k];
|
---|
| 260 | valP = new double[n * k];
|
---|
| 261 | var curP = new double[n - 1];
|
---|
| 262 | rowP[0] = 0;
|
---|
[14806] | 263 | for(var i = 0; i < n; i++) rowP[i + 1] = rowP[i] + k;
|
---|
[14512] | 264 |
|
---|
[14788] | 265 | var objX = new List<IndexedItem<T>>();
|
---|
[14806] | 266 | for(var i = 0; i < n; i++) objX.Add(new IndexedItem<T>(i, x[i]));
|
---|
[14512] | 267 |
|
---|
[14788] | 268 | // Build ball tree on data set
|
---|
| 269 | var tree = new VantagePointTree<IndexedItem<T>>(new IndexedItemDistance<T>(distance), objX); // do we really want to re-create the tree on each call?
|
---|
[14742] | 270 |
|
---|
[14788] | 271 | // Loop over all points to find nearest neighbors
|
---|
[14806] | 272 | for(var i = 0; i < n; i++) {
|
---|
[14788] | 273 | IList<IndexedItem<T>> indices;
|
---|
| 274 | IList<double> distances;
|
---|
[14742] | 275 |
|
---|
[14788] | 276 | // Find nearest neighbors
|
---|
| 277 | tree.Search(objX[i], k + 1, out indices, out distances);
|
---|
[14512] | 278 |
|
---|
[14788] | 279 | // Initialize some variables for binary search
|
---|
| 280 | var found = false;
|
---|
| 281 | var beta = 1.0;
|
---|
| 282 | var minBeta = double.MinValue;
|
---|
| 283 | var maxBeta = double.MaxValue;
|
---|
| 284 | const double tol = 1e-5; // TODO: why 1e-5?
|
---|
[14512] | 285 |
|
---|
[14788] | 286 | // Iterate until we found a good perplexity
|
---|
| 287 | var iter = 0; double sumP = 0;
|
---|
[14806] | 288 | while(!found && iter < 200) { // TODO 200 iterations always ok?
|
---|
[14512] | 289 |
|
---|
[14788] | 290 | // Compute Gaussian kernel row
|
---|
[14806] | 291 | for(var m = 0; m < k; m++) curP[m] = Math.Exp(-beta * distances[m + 1]); // TODO distances m+1?
|
---|
[14512] | 292 |
|
---|
[14788] | 293 | // Compute entropy of current row
|
---|
| 294 | sumP = double.Epsilon;
|
---|
[14806] | 295 | for(var m = 0; m < k; m++) sumP += curP[m];
|
---|
[14788] | 296 | var h = .0;
|
---|
[14806] | 297 | for(var m = 0; m < k; m++) h += beta * (distances[m + 1] * curP[m]); // TODO: distances m+1?
|
---|
[14788] | 298 | h = h / sumP + Math.Log(sumP);
|
---|
| 299 |
|
---|
| 300 | // Evaluate whether the entropy is within the tolerance level
|
---|
| 301 | var hdiff = h - Math.Log(perplexity);
|
---|
[14806] | 302 | if(hdiff < tol && -hdiff < tol) {
|
---|
[14788] | 303 | found = true;
|
---|
| 304 | } else {
|
---|
[14806] | 305 | if(hdiff > 0) {
|
---|
[14788] | 306 | minBeta = beta;
|
---|
[14806] | 307 | if(maxBeta.IsAlmost(double.MaxValue) || maxBeta.IsAlmost(double.MinValue))
|
---|
[14788] | 308 | beta *= 2.0;
|
---|
| 309 | else
|
---|
| 310 | beta = (beta + maxBeta) / 2.0;
|
---|
| 311 | } else {
|
---|
| 312 | maxBeta = beta;
|
---|
[14806] | 313 | if(minBeta.IsAlmost(double.MinValue) || minBeta.IsAlmost(double.MaxValue))
|
---|
[14788] | 314 | beta /= 2.0;
|
---|
| 315 | else
|
---|
| 316 | beta = (beta + minBeta) / 2.0;
|
---|
| 317 | }
|
---|
| 318 | }
|
---|
| 319 |
|
---|
| 320 | // Update iteration counter
|
---|
| 321 | iter++;
|
---|
| 322 | }
|
---|
| 323 |
|
---|
| 324 | // Row-normalize current row of P and store in matrix
|
---|
[14806] | 325 | for(var m = 0; m < k; m++) curP[m] /= sumP;
|
---|
| 326 | for(var m = 0; m < k; m++) {
|
---|
[14788] | 327 | colP[rowP[i] + m] = indices[m + 1].Index;
|
---|
| 328 | valP[rowP[i] + m] = curP[m];
|
---|
| 329 | }
|
---|
[14512] | 330 | }
|
---|
| 331 | }
|
---|
[14788] | 332 | private static void ComputeGaussianPerplexity(T[] x, IDistance<T> distance, double[,] p, double perplexity) {
|
---|
| 333 | // Compute the distance matrix
|
---|
| 334 | var dd = ComputeDistances(x, distance);
|
---|
| 335 |
|
---|
| 336 | int n = x.Length;
|
---|
| 337 | // Compute the Gaussian kernel row by row
|
---|
[14806] | 338 | for(var i = 0; i < n; i++) {
|
---|
[14788] | 339 | // Initialize some variables
|
---|
| 340 | var found = false;
|
---|
| 341 | var beta = 1.0;
|
---|
| 342 | var minBeta = -double.MaxValue;
|
---|
| 343 | var maxBeta = double.MaxValue;
|
---|
| 344 | const double tol = 1e-5;
|
---|
| 345 | double sumP = 0;
|
---|
| 346 |
|
---|
| 347 | // Iterate until we found a good perplexity
|
---|
| 348 | var iter = 0;
|
---|
[14806] | 349 | while(!found && iter < 200) { // TODO constant
|
---|
[14788] | 350 |
|
---|
| 351 | // Compute Gaussian kernel row
|
---|
[14806] | 352 | for(var m = 0; m < n; m++) p[i, m] = Math.Exp(-beta * dd[i][m]);
|
---|
[14788] | 353 | p[i, i] = double.Epsilon;
|
---|
| 354 |
|
---|
| 355 | // Compute entropy of current row
|
---|
| 356 | sumP = double.Epsilon;
|
---|
[14806] | 357 | for(var m = 0; m < n; m++) sumP += p[i, m];
|
---|
[14788] | 358 | var h = 0.0;
|
---|
[14806] | 359 | for(var m = 0; m < n; m++) h += beta * (dd[i][m] * p[i, m]);
|
---|
[14788] | 360 | h = h / sumP + Math.Log(sumP);
|
---|
| 361 |
|
---|
| 362 | // Evaluate whether the entropy is within the tolerance level
|
---|
| 363 | var hdiff = h - Math.Log(perplexity);
|
---|
[14806] | 364 | if(hdiff < tol && -hdiff < tol) {
|
---|
[14788] | 365 | found = true;
|
---|
| 366 | } else {
|
---|
[14806] | 367 | if(hdiff > 0) {
|
---|
[14788] | 368 | minBeta = beta;
|
---|
[14806] | 369 | if(maxBeta.IsAlmost(double.MaxValue) || maxBeta.IsAlmost(double.MinValue))
|
---|
[14788] | 370 | beta *= 2.0;
|
---|
| 371 | else
|
---|
| 372 | beta = (beta + maxBeta) / 2.0;
|
---|
| 373 | } else {
|
---|
| 374 | maxBeta = beta;
|
---|
[14806] | 375 | if(minBeta.IsAlmost(double.MinValue) || minBeta.IsAlmost(double.MaxValue))
|
---|
[14788] | 376 | beta /= 2.0;
|
---|
| 377 | else
|
---|
| 378 | beta = (beta + minBeta) / 2.0;
|
---|
| 379 | }
|
---|
| 380 | }
|
---|
| 381 |
|
---|
| 382 | // Update iteration counter
|
---|
| 383 | iter++;
|
---|
| 384 | }
|
---|
| 385 |
|
---|
| 386 | // Row normalize P
|
---|
[14806] | 387 | for(var m = 0; m < n; m++) p[i, m] /= sumP;
|
---|
[14512] | 388 | }
|
---|
| 389 | }
|
---|
| 390 |
|
---|
[14788] | 391 | private static double[][] ComputeDistances(T[] x, IDistance<T> distance) {
|
---|
[14806] | 392 | var res = new double[x.Length][];
|
---|
| 393 | for(int r = 0; r < x.Length; r++) {
|
---|
| 394 | var rowV = new double[x.Length];
|
---|
| 395 | // all distances must be symmetric
|
---|
| 396 | for(int c = 0; c < r; c++) {
|
---|
| 397 | rowV[c] = res[c][r];
|
---|
| 398 | }
|
---|
| 399 | rowV[r] = 0.0; // distance to self is zero for all distances
|
---|
| 400 | for(int c = r + 1; c < x.Length; c++) {
|
---|
| 401 | rowV[c] = distance.Get(x[r], x[c]);
|
---|
| 402 | }
|
---|
| 403 | res[r] = rowV;
|
---|
| 404 | }
|
---|
| 405 | return res;
|
---|
| 406 | // return x.Select(m => x.Select(n => distance.Get(m, n)).ToArray()).ToArray();
|
---|
[14788] | 407 | }
|
---|
[14414] | 408 |
|
---|
[14788] | 409 | private static double EvaluateErrorExact(double[,] p, double[,] y, int n, int d) {
|
---|
| 410 | // Compute the squared Euclidean distance matrix
|
---|
| 411 | var dd = new double[n, n];
|
---|
| 412 | var q = new double[n, n];
|
---|
[14806] | 413 | ComputeSquaredEuclideanDistance(y, n, d, dd); // TODO: we use Euclidian distance regardless of the actual distance function
|
---|
[14414] | 414 |
|
---|
[14788] | 415 | // Compute Q-matrix and normalization sum
|
---|
| 416 | var sumQ = double.Epsilon;
|
---|
[14806] | 417 | for(var n1 = 0; n1 < n; n1++) {
|
---|
| 418 | for(var m = 0; m < n; m++) {
|
---|
| 419 | if(n1 != m) {
|
---|
[14788] | 420 | q[n1, m] = 1 / (1 + dd[n1, m]);
|
---|
| 421 | sumQ += q[n1, m];
|
---|
| 422 | } else q[n1, m] = double.Epsilon;
|
---|
| 423 | }
|
---|
| 424 | }
|
---|
[14806] | 425 | for(var i = 0; i < n; i++) for(var j = 0; j < n; j++) q[i, j] /= sumQ;
|
---|
[14414] | 426 |
|
---|
[14788] | 427 | // Sum t-SNE error
|
---|
| 428 | var c = .0;
|
---|
[14806] | 429 | for(var i = 0; i < n; i++)
|
---|
| 430 | for(var j = 0; j < n; j++) {
|
---|
[14788] | 431 | c += p[i, j] * Math.Log((p[i, j] + float.Epsilon) / (q[i, j] + float.Epsilon));
|
---|
| 432 | }
|
---|
| 433 | return c;
|
---|
| 434 | }
|
---|
[14806] | 435 |
|
---|
| 436 | // TODO: there seems to be a bug in the error approximation.
|
---|
| 437 | // The mapping of the approximate tSNE looks good but the error curve never changes.
|
---|
[14788] | 438 | private static double EvaluateErrorApproximate(IReadOnlyList<int> rowP, IReadOnlyList<int> colP, IReadOnlyList<double> valP, double[,] y, double theta) {
|
---|
| 439 | // Get estimate of normalization term
|
---|
| 440 | var n = y.GetLength(0);
|
---|
| 441 | var d = y.GetLength(1);
|
---|
| 442 | var tree = new SpacePartitioningTree(y);
|
---|
| 443 | var buff = new double[d];
|
---|
| 444 | double sumQ = 0.0;
|
---|
[14806] | 445 | for(var i = 0; i < n; i++) tree.ComputeNonEdgeForces(i, theta, buff, ref sumQ);
|
---|
[14414] | 446 |
|
---|
[14788] | 447 | // Loop over all edges to compute t-SNE error
|
---|
| 448 | var c = .0;
|
---|
[14806] | 449 | for(var k = 0; k < n; k++) {
|
---|
| 450 | for(var i = rowP[k]; i < rowP[k + 1]; i++) {
|
---|
[14788] | 451 | var q = .0;
|
---|
[14806] | 452 | for(var j = 0; j < d; j++) buff[j] = y[k, j];
|
---|
| 453 | for(var j = 0; j < d; j++) buff[j] -= y[colP[i], j];
|
---|
| 454 | for(var j = 0; j < d; j++) q += buff[j] * buff[j]; // TODO: squared error is used here!
|
---|
[14788] | 455 | q = 1.0 / (1.0 + q) / sumQ;
|
---|
| 456 | c += valP[i] * Math.Log((valP[i] + float.Epsilon) / (q + float.Epsilon));
|
---|
| 457 | }
|
---|
| 458 | }
|
---|
| 459 | return c;
|
---|
| 460 | }
|
---|
| 461 | private static void SymmetrizeMatrix(IReadOnlyList<int> rowP, IReadOnlyList<int> colP, IReadOnlyList<double> valP, out int[] symRowP, out int[] symColP, out double[] symValP) {
|
---|
[14414] | 462 |
|
---|
[14788] | 463 | // Count number of elements and row counts of symmetric matrix
|
---|
| 464 | var n = rowP.Count - 1;
|
---|
| 465 | var rowCounts = new int[n];
|
---|
[14806] | 466 | for(var j = 0; j < n; j++) {
|
---|
| 467 | for(var i = rowP[j]; i < rowP[j + 1]; i++) {
|
---|
[14414] | 468 |
|
---|
[14788] | 469 | // Check whether element (col_P[i], n) is present
|
---|
| 470 | var present = false;
|
---|
[14806] | 471 | for(var m = rowP[colP[i]]; m < rowP[colP[i] + 1]; m++) {
|
---|
| 472 | if(colP[m] == j) present = true;
|
---|
[14788] | 473 | }
|
---|
[14806] | 474 | if(present) rowCounts[j]++;
|
---|
[14788] | 475 | else {
|
---|
| 476 | rowCounts[j]++;
|
---|
| 477 | rowCounts[colP[i]]++;
|
---|
| 478 | }
|
---|
| 479 | }
|
---|
| 480 | }
|
---|
| 481 | var noElem = 0;
|
---|
[14806] | 482 | for(var i = 0; i < n; i++) noElem += rowCounts[i];
|
---|
[14414] | 483 |
|
---|
[14788] | 484 | // Allocate memory for symmetrized matrix
|
---|
| 485 | symRowP = new int[n + 1];
|
---|
| 486 | symColP = new int[noElem];
|
---|
| 487 | symValP = new double[noElem];
|
---|
[14414] | 488 |
|
---|
[14788] | 489 | // Construct new row indices for symmetric matrix
|
---|
| 490 | symRowP[0] = 0;
|
---|
[14806] | 491 | for(var i = 0; i < n; i++) symRowP[i + 1] = symRowP[i] + rowCounts[i];
|
---|
[14788] | 492 |
|
---|
| 493 | // Fill the result matrix
|
---|
| 494 | var offset = new int[n];
|
---|
[14806] | 495 | for(var j = 0; j < n; j++) {
|
---|
| 496 | for(var i = rowP[j]; i < rowP[j + 1]; i++) { // considering element(n, colP[i])
|
---|
[14788] | 497 |
|
---|
| 498 | // Check whether element (col_P[i], n) is present
|
---|
| 499 | var present = false;
|
---|
[14806] | 500 | for(var m = rowP[colP[i]]; m < rowP[colP[i] + 1]; m++) {
|
---|
| 501 | if(colP[m] != j) continue;
|
---|
[14788] | 502 | present = true;
|
---|
[14806] | 503 | if(j > colP[i]) continue; // make sure we do not add elements twice
|
---|
[14788] | 504 | symColP[symRowP[j] + offset[j]] = colP[i];
|
---|
| 505 | symColP[symRowP[colP[i]] + offset[colP[i]]] = j;
|
---|
| 506 | symValP[symRowP[j] + offset[j]] = valP[i] + valP[m];
|
---|
| 507 | symValP[symRowP[colP[i]] + offset[colP[i]]] = valP[i] + valP[m];
|
---|
[14414] | 508 | }
|
---|
[14788] | 509 |
|
---|
| 510 | // If (colP[i], n) is not present, there is no addition involved
|
---|
[14806] | 511 | if(!present) {
|
---|
[14788] | 512 | symColP[symRowP[j] + offset[j]] = colP[i];
|
---|
| 513 | symColP[symRowP[colP[i]] + offset[colP[i]]] = j;
|
---|
| 514 | symValP[symRowP[j] + offset[j]] = valP[i];
|
---|
| 515 | symValP[symRowP[colP[i]] + offset[colP[i]]] = valP[i];
|
---|
| 516 | }
|
---|
| 517 |
|
---|
| 518 | // Update offsets
|
---|
[14806] | 519 | if(present && (j > colP[i])) continue;
|
---|
[14788] | 520 | offset[j]++;
|
---|
[14806] | 521 | if(colP[i] != j) offset[colP[i]]++;
|
---|
[14414] | 522 | }
|
---|
| 523 | }
|
---|
| 524 |
|
---|
[14806] | 525 | for(var i = 0; i < noElem; i++) symValP[i] /= 2.0;
|
---|
[14414] | 526 | }
|
---|
[14788] | 527 |
|
---|
[14414] | 528 | }
|
---|
[14785] | 529 |
|
---|
[14788] | 530 | public static TSNEState CreateState(T[] data, IDistance<T> distance, IRandom random, int newDimensions = 2, double perplexity = 25, double theta = 0,
|
---|
| 531 | int stopLyingIter = 250, int momSwitchIter = 250, double momentum = .5, double finalMomentum = .8, double eta = 200.0
|
---|
| 532 | ) {
|
---|
| 533 | return new TSNEState(data, distance, random, newDimensions, perplexity, theta, stopLyingIter, momSwitchIter, momentum, finalMomentum, eta);
|
---|
| 534 | }
|
---|
[14414] | 535 |
|
---|
| 536 |
|
---|
[14788] | 537 | public static double[,] Iterate(TSNEState state) {
|
---|
[14806] | 538 | if(state.exact)
|
---|
[14788] | 539 | ComputeExactGradient(state.p, state.newData, state.noDatapoints, state.newDimensions, state.dY);
|
---|
| 540 | else
|
---|
| 541 | ComputeApproximateGradient(state.rowP, state.colP, state.valP, state.newData, state.noDatapoints, state.newDimensions, state.dY, state.theta);
|
---|
[14414] | 542 |
|
---|
[14788] | 543 | // Update gains
|
---|
[14806] | 544 | for(var i = 0; i < state.noDatapoints; i++) {
|
---|
| 545 | for(var j = 0; j < state.newDimensions; j++) {
|
---|
[14788] | 546 | state.gains[i, j] = Math.Sign(state.dY[i, j]) != Math.Sign(state.uY[i, j])
|
---|
| 547 | ? state.gains[i, j] + .2
|
---|
| 548 | : state.gains[i, j] * .8; // 20% up or 20% down // TODO: +0.2?!
|
---|
[14414] | 549 |
|
---|
[14806] | 550 | if(state.gains[i, j] < .01) state.gains[i, j] = .01; // TODO why limit the gains?
|
---|
[14414] | 551 | }
|
---|
[14788] | 552 | }
|
---|
[14414] | 553 |
|
---|
[14788] | 554 |
|
---|
| 555 | // Perform gradient update (with momentum and gains)
|
---|
[14806] | 556 | for(var i = 0; i < state.noDatapoints; i++)
|
---|
| 557 | for(var j = 0; j < state.newDimensions; j++)
|
---|
[14788] | 558 | state.uY[i, j] = state.currentMomentum * state.uY[i, j] - state.eta * state.gains[i, j] * state.dY[i, j];
|
---|
| 559 |
|
---|
[14806] | 560 | for(var i = 0; i < state.noDatapoints; i++)
|
---|
| 561 | for(var j = 0; j < state.newDimensions; j++)
|
---|
[14788] | 562 | state.newData[i, j] = state.newData[i, j] + state.uY[i, j];
|
---|
| 563 |
|
---|
| 564 | // Make solution zero-mean
|
---|
| 565 | ZeroMean(state.newData);
|
---|
| 566 | // Stop lying about the P-values after a while, and switch momentum
|
---|
| 567 |
|
---|
[14806] | 568 | if(state.iter == state.stopLyingIter) {
|
---|
| 569 | if(state.exact)
|
---|
| 570 | 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?
|
---|
[14788] | 571 | else
|
---|
[14806] | 572 | for(var i = 0; i < state.rowP[state.noDatapoints]; i++) state.valP[i] /= 12.0; // XXX are we not scaling all values?
|
---|
[14414] | 573 | }
|
---|
[14788] | 574 |
|
---|
[14806] | 575 | if(state.iter == state.momSwitchIter)
|
---|
[14788] | 576 | state.currentMomentum = state.finalMomentum;
|
---|
| 577 |
|
---|
| 578 | state.iter++;
|
---|
| 579 | return state.newData;
|
---|
[14414] | 580 | }
|
---|
[14785] | 581 |
|
---|
[14788] | 582 |
|
---|
| 583 | private static void ComputeApproximateGradient(int[] rowP, int[] colP, double[] valP, double[,] y, int n, int d, double[,] dC, double theta) {
|
---|
| 584 | var tree = new SpacePartitioningTree(y);
|
---|
| 585 | double sumQ = 0.0;
|
---|
| 586 | var posF = new double[n, d];
|
---|
| 587 | var negF = new double[n, d];
|
---|
| 588 | tree.ComputeEdgeForces(rowP, colP, valP, n, posF);
|
---|
| 589 | var row = new double[d];
|
---|
[14806] | 590 | for(var n1 = 0; n1 < n; n1++) {
|
---|
[14788] | 591 | Buffer.BlockCopy(negF, (sizeof(double) * n1 * d), row, 0, d);
|
---|
| 592 | tree.ComputeNonEdgeForces(n1, theta, row, ref sumQ);
|
---|
| 593 | }
|
---|
| 594 |
|
---|
| 595 | // Compute final t-SNE gradient
|
---|
[14806] | 596 | for(var i = 0; i < n; i++)
|
---|
| 597 | for(var j = 0; j < d; j++) {
|
---|
[14788] | 598 | dC[i, j] = posF[i, j] - negF[i, j] / sumQ;
|
---|
| 599 | }
|
---|
[14414] | 600 | }
|
---|
[14785] | 601 |
|
---|
[14414] | 602 | private static void ComputeExactGradient(double[,] p, double[,] y, int n, int d, double[,] dC) {
|
---|
| 603 |
|
---|
| 604 | // Make sure the current gradient contains zeros
|
---|
[14806] | 605 | for(var i = 0; i < n; i++) for(var j = 0; j < d; j++) dC[i, j] = 0.0;
|
---|
[14414] | 606 |
|
---|
| 607 | // Compute the squared Euclidean distance matrix
|
---|
| 608 | var dd = new double[n, n];
|
---|
[14806] | 609 | ComputeSquaredEuclideanDistance(y, n, d, dd); // TODO: we use Euclidian distance regardless which distance function is actually set!
|
---|
[14414] | 610 |
|
---|
| 611 | // Compute Q-matrix and normalization sum
|
---|
| 612 | var q = new double[n, n];
|
---|
| 613 | var sumQ = .0;
|
---|
[14806] | 614 | for(var n1 = 0; n1 < n; n1++) {
|
---|
| 615 | for(var m = 0; m < n; m++) {
|
---|
| 616 | if(n1 == m) continue;
|
---|
[14414] | 617 | q[n1, m] = 1 / (1 + dd[n1, m]);
|
---|
| 618 | sumQ += q[n1, m];
|
---|
| 619 | }
|
---|
| 620 | }
|
---|
| 621 |
|
---|
| 622 | // Perform the computation of the gradient
|
---|
[14806] | 623 | for(var n1 = 0; n1 < n; n1++) {
|
---|
| 624 | for(var m = 0; m < n; m++) {
|
---|
| 625 | if(n1 == m) continue;
|
---|
[14414] | 626 | var mult = (p[n1, m] - q[n1, m] / sumQ) * q[n1, m];
|
---|
[14806] | 627 | for(var d1 = 0; d1 < d; d1++) {
|
---|
[14414] | 628 | dC[n1, d1] += (y[n1, d1] - y[m, d1]) * mult;
|
---|
| 629 | }
|
---|
| 630 | }
|
---|
| 631 | }
|
---|
| 632 | }
|
---|
[14788] | 633 |
|
---|
[14414] | 634 | private static void ComputeSquaredEuclideanDistance(double[,] x, int n, int d, double[,] dd) {
|
---|
| 635 | var dataSums = new double[n];
|
---|
[14806] | 636 | for(var i = 0; i < n; i++) {
|
---|
| 637 | for(var j = 0; j < d; j++) {
|
---|
[14414] | 638 | dataSums[i] += x[i, j] * x[i, j];
|
---|
| 639 | }
|
---|
| 640 | }
|
---|
[14806] | 641 | for(var i = 0; i < n; i++) {
|
---|
| 642 | for(var m = 0; m < n; m++) {
|
---|
[14414] | 643 | dd[i, m] = dataSums[i] + dataSums[m];
|
---|
| 644 | }
|
---|
| 645 | }
|
---|
[14806] | 646 | for(var i = 0; i < n; i++) {
|
---|
[14414] | 647 | dd[i, i] = 0.0;
|
---|
[14806] | 648 | for(var m = i + 1; m < n; m++) {
|
---|
[14414] | 649 | dd[i, m] = 0.0;
|
---|
[14806] | 650 | for(var j = 0; j < d; j++) {
|
---|
[14414] | 651 | dd[i, m] += (x[i, j] - x[m, j]) * (x[i, j] - x[m, j]);
|
---|
| 652 | }
|
---|
| 653 | dd[m, i] = dd[i, m];
|
---|
| 654 | }
|
---|
| 655 | }
|
---|
| 656 | }
|
---|
| 657 |
|
---|
| 658 | private static void ZeroMean(double[,] x) {
|
---|
| 659 | // Compute data mean
|
---|
| 660 | var n = x.GetLength(0);
|
---|
| 661 | var d = x.GetLength(1);
|
---|
| 662 | var mean = new double[d];
|
---|
[14806] | 663 | for(var i = 0; i < n; i++) {
|
---|
| 664 | for(var j = 0; j < d; j++) {
|
---|
[14414] | 665 | mean[j] += x[i, j];
|
---|
| 666 | }
|
---|
| 667 | }
|
---|
[14806] | 668 | for(var i = 0; i < d; i++) {
|
---|
[14414] | 669 | mean[i] /= n;
|
---|
| 670 | }
|
---|
| 671 | // Subtract data mean
|
---|
[14806] | 672 | for(var i = 0; i < n; i++) {
|
---|
| 673 | for(var j = 0; j < d; j++) {
|
---|
[14414] | 674 | x[i, j] -= mean[j];
|
---|
| 675 | }
|
---|
| 676 | }
|
---|
| 677 | }
|
---|
| 678 | }
|
---|
| 679 | }
|
---|