Free cookie consent management tool by TermsFeed Policy Generator

source: stable/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/TSNEStatic.cs @ 17529

Last change on this file since 17529 was 17181, checked in by swagner, 5 years ago

#2875: Merged r17180 from trunk to stable

File size: 28.9 KB
Line 
1#region License Information
2/* HeuristicLab
3 * Copyright (C) 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
56using System;
57using System.Collections.Generic;
58using HeuristicLab.Collections;
59using HeuristicLab.Common;
60using HeuristicLab.Core;
61using HEAL.Attic;
62using HeuristicLab.Random;
63
64namespace HeuristicLab.Algorithms.DataAnalysis {
65  [StorableType("23E56F69-7AA1-4853-BFEC-7B4DBC346B47")]
66  public class TSNEStatic<T> {
67    [StorableConstructor]
68    protected TSNEStatic(StorableConstructorFlag _) {
69    }
70
71    [StorableType("1878DFC2-5DDF-4BBD-85BA-4A6E19AB59C2")]
72    public sealed class TSNEState : DeepCloneable {
73      #region Storables
74      // initialized once
75      [Storable]
76      public IDistance<T> distance;
77      [Storable]
78      public IRandom random;
79      [Storable]
80      public double perplexity;
81      [Storable]
82      public bool exact;
83      [Storable]
84      public int noDatapoints;
85      [Storable]
86      public double finalMomentum;
87      [Storable]
88      public int momSwitchIter;
89      [Storable]
90      public int stopLyingIter;
91      [Storable]
92      public double theta;
93      [Storable]
94      public double eta;
95      [Storable]
96      public int newDimensions;
97
98      // for approximate version: sparse representation of similarity/distance matrix
99      [Storable]
100      public double[] valP; // similarity/distance
101      [Storable]
102      public int[] rowP; // row index
103      [Storable]
104      public int[] colP; // col index
105
106      // for exact version: dense representation of distance/similarity matrix
107      [Storable]
108      public double[,] p;
109
110      // mapped data
111      [Storable]
112      public double[,] newData;
113
114      [Storable]
115      public int iter;
116      [Storable]
117      public double currentMomentum;
118
119      // helper variables (updated in each iteration)
120      [Storable]
121      public double[,] gains;
122      [Storable]
123      public double[,] uY;
124      [Storable]
125      public double[,] dY;
126      #endregion
127
128      #region Constructors & Cloning
129      private TSNEState(TSNEState original, Cloner cloner) : base(original, cloner) {
130        distance = cloner.Clone(original.distance);
131        random = cloner.Clone(original.random);
132        perplexity = original.perplexity;
133        exact = original.exact;
134        noDatapoints = original.noDatapoints;
135        finalMomentum = original.finalMomentum;
136        momSwitchIter = original.momSwitchIter;
137        stopLyingIter = original.stopLyingIter;
138        theta = original.theta;
139        eta = original.eta;
140        newDimensions = original.newDimensions;
141        if (original.valP != null) {
142          valP = new double[original.valP.Length];
143          Array.Copy(original.valP, valP, valP.Length);
144        }
145        if (original.rowP != null) {
146          rowP = new int[original.rowP.Length];
147          Array.Copy(original.rowP, rowP, rowP.Length);
148        }
149        if (original.colP != null) {
150          colP = new int[original.colP.Length];
151          Array.Copy(original.colP, colP, colP.Length);
152        }
153        if (original.p != null) {
154          p = new double[original.p.GetLength(0), original.p.GetLength(1)];
155          Array.Copy(original.p, p, p.Length);
156        }
157        newData = new double[original.newData.GetLength(0), original.newData.GetLength(1)];
158        Array.Copy(original.newData, newData, newData.Length);
159        iter = original.iter;
160        currentMomentum = original.currentMomentum;
161        gains = new double[original.gains.GetLength(0), original.gains.GetLength(1)];
162        Array.Copy(original.gains, gains, gains.Length);
163        uY = new double[original.uY.GetLength(0), original.uY.GetLength(1)];
164        Array.Copy(original.uY, uY, uY.Length);
165        dY = new double[original.dY.GetLength(0), original.dY.GetLength(1)];
166        Array.Copy(original.dY, dY, dY.Length);
167      }
168
169      public override IDeepCloneable Clone(Cloner cloner) {
170        return new TSNEState(this, cloner);
171      }
172
173      [StorableConstructor]
174      private TSNEState(StorableConstructorFlag _) { }
175
176      public TSNEState(IReadOnlyList<T> data, IDistance<T> distance, IRandom random, int newDimensions, double perplexity,
177        double theta, int stopLyingIter, int momSwitchIter, double momentum, double finalMomentum, double eta, bool randomInit) {
178        this.distance = distance;
179        this.random = random;
180        this.newDimensions = newDimensions;
181        this.perplexity = perplexity;
182        this.theta = theta;
183        this.stopLyingIter = stopLyingIter;
184        this.momSwitchIter = momSwitchIter;
185        currentMomentum = momentum;
186        this.finalMomentum = finalMomentum;
187        this.eta = eta;
188
189        // initialize
190        noDatapoints = data.Count;
191        if (noDatapoints - 1 < 3 * perplexity)
192          throw new ArgumentException("Perplexity too large for the number of data points!");
193
194        exact = Math.Abs(theta) < double.Epsilon;
195        newData = new double[noDatapoints, newDimensions];
196        dY = new double[noDatapoints, newDimensions];
197        uY = new double[noDatapoints, newDimensions];
198        gains = new double[noDatapoints, newDimensions];
199        for (var i = 0; i < noDatapoints; i++)
200        for (var j = 0; j < newDimensions; j++)
201          gains[i, j] = 1.0;
202
203        p = null;
204        rowP = null;
205        colP = null;
206        valP = null;
207
208        //Calculate Similarities
209        if (exact) p = CalculateExactSimilarites(data, distance, perplexity);
210        else CalculateApproximateSimilarities(data, distance, perplexity, out rowP, out colP, out valP);
211
212        // Lie about the P-values (factor is 4 in the MATLAB implementation)
213        if (exact) for (var i = 0; i < noDatapoints; i++) for (var j = 0; j < noDatapoints; j++) p[i, j] *= 12.0;
214        else for (var i = 0; i < rowP[noDatapoints]; i++) valP[i] *= 12.0;
215
216        // Initialize solution (randomly)
217        var rand = new NormalDistributedRandom(random, 0, 1);
218        for (var i = 0; i < noDatapoints; i++)
219        for (var j = 0; j < newDimensions; j++)
220          newData[i, j] = rand.NextDouble() * .0001;
221
222        if (!(data[0] is IReadOnlyList<double>) || randomInit) return;
223        for (var i = 0; i < noDatapoints; i++)
224        for (var j = 0; j < newDimensions; j++) {
225          var row = (IReadOnlyList<double>) data[i];
226          newData[i, j] = row[j % row.Count];
227        }
228      }
229      #endregion
230
231      public double EvaluateError() {
232        return exact ? EvaluateErrorExact(p, newData, noDatapoints, newDimensions) : EvaluateErrorApproximate(rowP, colP, valP, newData, theta);
233      }
234
235      #region Helpers
236      private static void CalculateApproximateSimilarities(IReadOnlyList<T> data, IDistance<T> distance, double perplexity, out int[] rowP, out int[] colP, out double[] valP) {
237        // Compute asymmetric pairwise input similarities
238        ComputeGaussianPerplexity(data, distance, out rowP, out colP, out valP, perplexity, (int) (3 * perplexity));
239        // Symmetrize input similarities
240        int[] sRowP, symColP;
241        double[] sValP;
242        SymmetrizeMatrix(rowP, colP, valP, out sRowP, out symColP, out sValP);
243        rowP = sRowP;
244        colP = symColP;
245        valP = sValP;
246        var sumP = .0;
247        for (var i = 0; i < rowP[data.Count]; i++) sumP += valP[i];
248        for (var i = 0; i < rowP[data.Count]; i++) valP[i] /= sumP;
249      }
250      private static double[,] CalculateExactSimilarites(IReadOnlyList<T> data, IDistance<T> distance, double perplexity) {
251        // Compute similarities
252        var p = new double[data.Count, data.Count];
253        ComputeGaussianPerplexity(data, distance, p, perplexity);
254        // Symmetrize input similarities
255        for (var n = 0; n < data.Count; n++) {
256          for (var m = n + 1; m < data.Count; m++) {
257            p[n, m] += p[m, n];
258            p[m, n] = p[n, m];
259          }
260        }
261        var sumP = .0;
262        for (var i = 0; i < data.Count; i++) {
263          for (var j = 0; j < data.Count; j++) {
264            sumP += p[i, j];
265          }
266        }
267        for (var i = 0; i < data.Count; i++) {
268          for (var j = 0; j < data.Count; j++) {
269            p[i, j] /= sumP;
270          }
271        }
272        return p;
273      }
274      private static void ComputeGaussianPerplexity(IReadOnlyList<T> x, IDistance<T> distance, out int[] rowP, out int[] colP, out double[] valP, double perplexity, int k) {
275        if (perplexity > k) throw new ArgumentException("Perplexity should be lower than k!");
276
277        var n = x.Count;
278        // Allocate the memory we need
279        rowP = new int[n + 1];
280        colP = new int[n * k];
281        valP = new double[n * k];
282        var curP = new double[n - 1];
283        rowP[0] = 0;
284        for (var i = 0; i < n; i++) rowP[i + 1] = rowP[i] + k;
285
286        var objX = new List<IndexedItem<T>>();
287        for (var i = 0; i < n; i++) objX.Add(new IndexedItem<T>(i, x[i]));
288
289        // Build ball tree on data set
290        var tree = new VantagePointTree<IndexedItem<T>>(new IndexedItemDistance<T>(distance), objX);
291
292        // Loop over all points to find nearest neighbors
293        for (var i = 0; i < n; i++) {
294          IList<IndexedItem<T>> indices;
295          IList<double> distances;
296
297          // Find nearest neighbors
298          tree.Search(objX[i], k + 1, out indices, out distances);
299
300          // Initialize some variables for binary search
301          var found = false;
302          var beta = 1.0;
303          var minBeta = double.MinValue;
304          var maxBeta = double.MaxValue;
305          const double tol = 1e-5;
306
307          // Iterate until we found a good perplexity
308          var iter = 0;
309          double sumP = 0;
310          while (!found && iter < 200) {
311            // Compute Gaussian kernel row
312            for (var m = 0; m < k; m++) curP[m] = Math.Exp(-beta * distances[m + 1]);
313
314            // Compute entropy of current row
315            sumP = double.Epsilon;
316            for (var m = 0; m < k; m++) sumP += curP[m];
317            var h = .0;
318            for (var m = 0; m < k; m++) h += beta * (distances[m + 1] * curP[m]);
319            h = h / sumP + Math.Log(sumP);
320
321            // Evaluate whether the entropy is within the tolerance level
322            var hdiff = h - Math.Log(perplexity);
323            if (hdiff < tol && -hdiff < tol) {
324              found = true;
325            }
326            else {
327              if (hdiff > 0) {
328                minBeta = beta;
329                if (maxBeta.IsAlmost(double.MaxValue) || maxBeta.IsAlmost(double.MinValue))
330                  beta *= 2.0;
331                else
332                  beta = (beta + maxBeta) / 2.0;
333              }
334              else {
335                maxBeta = beta;
336                if (minBeta.IsAlmost(double.MinValue) || minBeta.IsAlmost(double.MaxValue))
337                  beta /= 2.0;
338                else
339                  beta = (beta + minBeta) / 2.0;
340              }
341            }
342
343            // Update iteration counter
344            iter++;
345          }
346
347          // Row-normalize current row of P and store in matrix
348          for (var m = 0; m < k; m++) curP[m] /= sumP;
349          for (var m = 0; m < k; m++) {
350            colP[rowP[i] + m] = indices[m + 1].Index;
351            valP[rowP[i] + m] = curP[m];
352          }
353        }
354      }
355      private static void ComputeGaussianPerplexity(IReadOnlyList<T> x, IDistance<T> distance, double[,] p, double perplexity) {
356        // Compute the distance matrix
357        var dd = ComputeDistances(x, distance);
358
359        var n = x.Count;
360        // Compute the Gaussian kernel row by row
361        for (var i = 0; i < n; i++) {
362          // Initialize some variables
363          var found = false;
364          var beta = 1.0;
365          var minBeta = double.MinValue;
366          var maxBeta = double.MaxValue;
367          const double tol = 1e-5;
368          double sumP = 0;
369
370          // Iterate until we found a good perplexity
371          var iter = 0;
372          while (!found && iter < 200) { // 200 iterations as in tSNE implementation by van der Maarten
373
374            // Compute Gaussian kernel row
375            for (var m = 0; m < n; m++) p[i, m] = Math.Exp(-beta * dd[i][m]);
376            p[i, i] = double.Epsilon;
377
378            // Compute entropy of current row
379            sumP = double.Epsilon;
380            for (var m = 0; m < n; m++) sumP += p[i, m];
381            var h = 0.0;
382            for (var m = 0; m < n; m++) h += beta * (dd[i][m] * p[i, m]);
383            h = h / sumP + Math.Log(sumP);
384
385            // Evaluate whether the entropy is within the tolerance level
386            var hdiff = h - Math.Log(perplexity);
387            if (hdiff < tol && -hdiff < tol) {
388              found = true;
389            }
390            else {
391              if (hdiff > 0) {
392                minBeta = beta;
393                if (maxBeta.IsAlmost(double.MaxValue) || maxBeta.IsAlmost(double.MinValue))
394                  beta *= 2.0;
395                else
396                  beta = (beta + maxBeta) / 2.0;
397              }
398              else {
399                maxBeta = beta;
400                if (minBeta.IsAlmost(double.MinValue) || minBeta.IsAlmost(double.MaxValue))
401                  beta /= 2.0;
402                else
403                  beta = (beta + minBeta) / 2.0;
404              }
405            }
406
407            // Update iteration counter
408            iter++;
409          }
410
411          // Row normalize P
412          for (var m = 0; m < n; m++) p[i, m] /= sumP;
413        }
414      }
415      private static double[][] ComputeDistances(IReadOnlyList<T> x, IDistance<T> distance) {
416        var res = new double[x.Count][];
417        for (var r = 0; r < x.Count; r++) {
418          var rowV = new double[x.Count];
419          // all distances must be symmetric
420          for (var c = 0; c < r; c++) {
421            rowV[c] = res[c][r];
422          }
423          rowV[r] = 0.0; // distance to self is zero for all distances
424          for (var c = r + 1; c < x.Count; c++) {
425            rowV[c] = distance.Get(x[r], x[c]);
426          }
427          res[r] = rowV;
428        }
429        return res;
430        // return x.Select(m => x.Select(n => distance.Get(m, n)).ToArray()).ToArray();
431      }
432      private static double EvaluateErrorExact(double[,] p, double[,] y, int n, int d) {
433        // Compute the squared Euclidean distance matrix
434        var dd = new double[n, n];
435        var q = new double[n, n];
436        ComputeSquaredEuclideanDistance(y, n, d, dd);
437
438        // Compute Q-matrix and normalization sum
439        var sumQ = double.Epsilon;
440        for (var n1 = 0; n1 < n; n1++) {
441          for (var m = 0; m < n; m++) {
442            if (n1 != m) {
443              q[n1, m] = 1 / (1 + dd[n1, m]);
444              sumQ += q[n1, m];
445            }
446            else q[n1, m] = double.Epsilon;
447          }
448        }
449        for (var i = 0; i < n; i++) for (var j = 0; j < n; j++) q[i, j] /= sumQ;
450
451        // Sum t-SNE error
452        var c = .0;
453        for (var i = 0; i < n; i++)
454        for (var j = 0; j < n; j++) {
455          c += p[i, j] * Math.Log((p[i, j] + float.Epsilon) / (q[i, j] + float.Epsilon));
456        }
457        return c;
458      }
459      private static double EvaluateErrorApproximate(IReadOnlyList<int> rowP, IReadOnlyList<int> colP, IReadOnlyList<double> valP, double[,] y, double theta) {
460        // Get estimate of normalization term
461        var n = y.GetLength(0);
462        var d = y.GetLength(1);
463        var tree = new SpacePartitioningTree(y);
464        var buff = new double[d];
465        var sumQ = 0.0;
466        for (var i = 0; i < n; i++) tree.ComputeNonEdgeForces(i, theta, buff, ref sumQ);
467
468        // Loop over all edges to compute t-SNE error
469        var c = .0;
470        for (var k = 0; k < n; k++) {
471          for (var i = rowP[k]; i < rowP[k + 1]; i++) {
472            var q = .0;
473            for (var j = 0; j < d; j++) buff[j] = y[k, j];
474            for (var j = 0; j < d; j++) buff[j] -= y[colP[i], j];
475            for (var j = 0; j < d; j++) q += buff[j] * buff[j];
476            q = (1.0 / (1.0 + q)) / sumQ;
477            c += valP[i] * Math.Log((valP[i] + float.Epsilon) / (q + float.Epsilon));
478          }
479        }
480        return c;
481      }
482      private static void SymmetrizeMatrix(IReadOnlyList<int> rowP, IReadOnlyList<int> colP, IReadOnlyList<double> valP, out int[] symRowP, out int[] symColP, out double[] symValP) {
483        // Count number of elements and row counts of symmetric matrix
484        var n = rowP.Count - 1;
485        var rowCounts = new int[n];
486        for (var j = 0; j < n; j++) {
487          for (var i = rowP[j]; i < rowP[j + 1]; i++) {
488            // Check whether element (col_P[i], n) is present
489            var present = false;
490            for (var m = rowP[colP[i]]; m < rowP[colP[i] + 1]; m++) {
491              if (colP[m] == j) present = true;
492            }
493            if (present) rowCounts[j]++;
494            else {
495              rowCounts[j]++;
496              rowCounts[colP[i]]++;
497            }
498          }
499        }
500        var noElem = 0;
501        for (var i = 0; i < n; i++) noElem += rowCounts[i];
502
503        // Allocate memory for symmetrized matrix
504        symRowP = new int[n + 1];
505        symColP = new int[noElem];
506        symValP = new double[noElem];
507
508        // Construct new row indices for symmetric matrix
509        symRowP[0] = 0;
510        for (var i = 0; i < n; i++) symRowP[i + 1] = symRowP[i] + rowCounts[i];
511
512        // Fill the result matrix
513        var offset = new int[n];
514        for (var j = 0; j < n; j++) {
515          for (var i = rowP[j]; i < rowP[j + 1]; i++) { // considering element(n, colP[i])
516
517            // Check whether element (col_P[i], n) is present
518            var present = false;
519            for (var m = rowP[colP[i]]; m < rowP[colP[i] + 1]; m++) {
520              if (colP[m] != j) continue;
521              present = true;
522              if (j > colP[i]) continue; // make sure we do not add elements twice
523              symColP[symRowP[j] + offset[j]] = colP[i];
524              symColP[symRowP[colP[i]] + offset[colP[i]]] = j;
525              symValP[symRowP[j] + offset[j]] = valP[i] + valP[m];
526              symValP[symRowP[colP[i]] + offset[colP[i]]] = valP[i] + valP[m];
527            }
528
529            // If (colP[i], n) is not present, there is no addition involved
530            if (!present) {
531              symColP[symRowP[j] + offset[j]] = colP[i];
532              symColP[symRowP[colP[i]] + offset[colP[i]]] = j;
533              symValP[symRowP[j] + offset[j]] = valP[i];
534              symValP[symRowP[colP[i]] + offset[colP[i]]] = valP[i];
535            }
536
537            // Update offsets
538            if (present && (j > colP[i])) continue;
539            offset[j]++;
540            if (colP[i] != j) offset[colP[i]]++;
541          }
542        }
543
544        for (var i = 0; i < noElem; i++) symValP[i] /= 2.0;
545      }
546      #endregion
547    }
548
549    /// <summary>
550    /// Static interface to tSNE
551    /// </summary>
552    /// <param name="data"></param>
553    /// <param name="distance">The distance function used to differentiate similar from non-similar points, e.g. Euclidean distance.</param>
554    /// <param name="random">Random number generator</param>
555    /// <param name="newDimensions">Dimensionality of projected space (usually 2 for easy visual analysis).</param>
556    /// <param name="perplexity">Perplexity parameter of tSNE. Comparable to k in a k-nearest neighbour algorithm. Recommended value is floor(number of points /3) or lower</param>
557    /// <param name="iterations">Maximum number of iterations for gradient descent.</param>
558    /// <param name="theta">Value describing how much appoximated gradients my differ from exact gradients. Set to 0 for exact calculation and in [0,1] otherwise. CAUTION: exact calculation of forces requires building a non-sparse N*N matrix where N is the number of data points. This may exceed memory limitations.</param>
559    /// <param name="stopLyingIter">Number of iterations after which p is no longer approximated.</param>
560    /// <param name="momSwitchIter">Number of iterations after which the momentum in the gradient descent is switched.</param>
561    /// <param name="momentum">The initial momentum in the gradient descent.</param>
562    /// <param name="finalMomentum">The final momentum in gradient descent (after momentum switch).</param>
563    /// <param name="eta">Gradient descent learning rate.</param>
564    /// <returns></returns>
565    public static double[,] Run(T[] data, IDistance<T> distance, IRandom random,
566      int newDimensions = 2, double perplexity = 25, int iterations = 1000,
567      double theta = 0, int stopLyingIter = 0, int momSwitchIter = 0, double momentum = .5,
568      double finalMomentum = .8, double eta = 10.0
569    ) {
570      var state = CreateState(data, distance, random, newDimensions, perplexity,
571        theta, stopLyingIter, momSwitchIter, momentum, finalMomentum, eta);
572
573      for (var i = 0; i < iterations - 1; i++) {
574        Iterate(state);
575      }
576      return Iterate(state);
577    }
578
579    public static TSNEState CreateState(T[] data, IDistance<T> distance, IRandom random,
580      int newDimensions = 2, double perplexity = 25, double theta = 0,
581      int stopLyingIter = 0, int momSwitchIter = 0, double momentum = .5,
582      double finalMomentum = .8, double eta = 10.0, bool randomInit = true
583    ) {
584      return new TSNEState(data, distance, random, newDimensions, perplexity, theta, stopLyingIter, momSwitchIter, momentum, finalMomentum, eta, randomInit);
585    }
586
587    public static double[,] Iterate(TSNEState state) {
588      if (state.exact)
589        ComputeExactGradient(state.p, state.newData, state.noDatapoints, state.newDimensions, state.dY);
590      else
591        ComputeApproximateGradient(state.rowP, state.colP, state.valP, state.newData, state.noDatapoints, state.newDimensions, state.dY, state.theta);
592
593      // Update gains
594      for (var i = 0; i < state.noDatapoints; i++) {
595        for (var j = 0; j < state.newDimensions; j++) {
596          state.gains[i, j] = Math.Sign(state.dY[i, j]) != Math.Sign(state.uY[i, j])
597            ? state.gains[i, j] + .2 // +0.2 nd *0.8 are used in two separate implementations of tSNE -> seems to be correct
598            : state.gains[i, j] * .8;
599          if (state.gains[i, j] < .01) state.gains[i, j] = .01;
600        }
601      }
602
603      // Perform gradient update (with momentum and gains)
604      for (var i = 0; i < state.noDatapoints; i++)
605      for (var j = 0; j < state.newDimensions; j++)
606        state.uY[i, j] = state.currentMomentum * state.uY[i, j] - state.eta * state.gains[i, j] * state.dY[i, j];
607
608      for (var i = 0; i < state.noDatapoints; i++)
609      for (var j = 0; j < state.newDimensions; j++)
610        state.newData[i, j] = state.newData[i, j] + state.uY[i, j];
611
612      // Make solution zero-mean
613      ZeroMean(state.newData);
614
615      // Stop lying about the P-values after a while, and switch momentum
616      if (state.iter == state.stopLyingIter) {
617        if (state.exact)
618          for (var i = 0; i < state.noDatapoints; i++)
619          for (var j = 0; j < state.noDatapoints; j++)
620            state.p[i, j] /= 12.0;
621        else
622          for (var i = 0; i < state.rowP[state.noDatapoints]; i++)
623            state.valP[i] /= 12.0;
624      }
625
626      if (state.iter == state.momSwitchIter)
627        state.currentMomentum = state.finalMomentum;
628
629      state.iter++;
630      return state.newData;
631    }
632
633    #region Helpers
634    private static void ComputeApproximateGradient(int[] rowP, int[] colP, double[] valP, double[,] y, int n, int d, double[,] dC, double theta) {
635      var tree = new SpacePartitioningTree(y);
636      var sumQ = 0.0;
637      var posF = new double[n, d];
638      var negF = new double[n, d];
639      SpacePartitioningTree.ComputeEdgeForces(rowP, colP, valP, n, posF, y, d);
640      var row = new double[d];
641      for (var n1 = 0; n1 < n; n1++) {
642        Array.Clear(row, 0, row.Length);
643        tree.ComputeNonEdgeForces(n1, theta, row, ref sumQ);
644        Buffer.BlockCopy(row, 0, negF, (sizeof(double) * n1 * d), d * sizeof(double));
645      }
646
647      // Compute final t-SNE gradient
648      for (var i = 0; i < n; i++)
649      for (var j = 0; j < d; j++) {
650        dC[i, j] = posF[i, j] - negF[i, j] / sumQ;
651      }
652    }
653
654    private static void ComputeExactGradient(double[,] p, double[,] y, int n, int d, double[,] dC) {
655      // Make sure the current gradient contains zeros
656      for (var i = 0; i < n; i++) for (var j = 0; j < d; j++) dC[i, j] = 0.0;
657
658      // Compute the squared Euclidean distance matrix
659      var dd = new double[n, n];
660      ComputeSquaredEuclideanDistance(y, n, d, dd);
661
662      // Compute Q-matrix and normalization sum
663      var q = new double[n, n];
664      var sumQ = .0;
665      for (var n1 = 0; n1 < n; n1++) {
666        for (var m = 0; m < n; m++) {
667          if (n1 == m) continue;
668          q[n1, m] = 1 / (1 + dd[n1, m]);
669          sumQ += q[n1, m];
670        }
671      }
672
673      // Perform the computation of the gradient
674      for (var n1 = 0; n1 < n; n1++) {
675        for (var m = 0; m < n; m++) {
676          if (n1 == m) continue;
677          var mult = (p[n1, m] - q[n1, m] / sumQ) * q[n1, m];
678          for (var d1 = 0; d1 < d; d1++) {
679            dC[n1, d1] += (y[n1, d1] - y[m, d1]) * mult;
680          }
681        }
682      }
683    }
684
685    private static void ComputeSquaredEuclideanDistance(double[,] x, int n, int d, double[,] dd) {
686      var dataSums = new double[n];
687      for (var i = 0; i < n; i++) {
688        for (var j = 0; j < d; j++) {
689          dataSums[i] += x[i, j] * x[i, j];
690        }
691      }
692      for (var i = 0; i < n; i++) {
693        for (var m = 0; m < n; m++) {
694          dd[i, m] = dataSums[i] + dataSums[m];
695        }
696      }
697      for (var i = 0; i < n; i++) {
698        dd[i, i] = 0.0;
699        for (var m = i + 1; m < n; m++) {
700          dd[i, m] = 0.0;
701          for (var j = 0; j < d; j++) {
702            dd[i, m] += (x[i, j] - x[m, j]) * (x[i, j] - x[m, j]);
703          }
704          dd[m, i] = dd[i, m];
705        }
706      }
707    }
708
709    private static void ZeroMean(double[,] x) {
710      // Compute data mean
711      var n = x.GetLength(0);
712      var d = x.GetLength(1);
713      var mean = new double[d];
714      for (var i = 0; i < n; i++) {
715        for (var j = 0; j < d; j++) {
716          mean[j] += x[i, j];
717        }
718      }
719      for (var i = 0; i < d; i++) {
720        mean[i] /= n;
721      }
722      // Subtract data mean
723      for (var i = 0; i < n; i++) {
724        for (var j = 0; j < d; j++) {
725          x[i, j] -= mean[j];
726        }
727      }
728    }
729    #endregion
730  }
731}
Note: See TracBrowser for help on using the repository browser.