Free cookie consent management tool by TermsFeed Policy Generator

source: trunk/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/TSNEStatic.cs @ 17313

Last change on this file since 17313 was 17180, checked in by swagner, 5 years ago

#2875: Removed years in copyrights

File size: 28.9 KB
RevLine 
[14414]1#region License Information
2/* HeuristicLab
[17180]3 * Copyright (C) Heuristic and Evolutionary Algorithms Laboratory (HEAL)
[14414]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;
[14785]58using HeuristicLab.Collections;
[14414]59using HeuristicLab.Common;
60using HeuristicLab.Core;
[16565]61using HEAL.Attic;
[14414]62using HeuristicLab.Random;
63
64namespace HeuristicLab.Algorithms.DataAnalysis {
[16565]65  [StorableType("23E56F69-7AA1-4853-BFEC-7B4DBC346B47")]
[14807]66  public class TSNEStatic<T> {
[16565]67    [StorableConstructor]
68    protected TSNEStatic(StorableConstructorFlag _) {
69    }
70
71    [StorableType("1878DFC2-5DDF-4BBD-85BA-4A6E19AB59C2")]
[14788]72    public sealed class TSNEState : DeepCloneable {
[15207]73      #region Storables
[14788]74      // initialized once
[14806]75      [Storable]
[14788]76      public IDistance<T> distance;
[14806]77      [Storable]
[14788]78      public IRandom random;
[14806]79      [Storable]
[14788]80      public double perplexity;
[14806]81      [Storable]
[14788]82      public bool exact;
[14806]83      [Storable]
[14788]84      public int noDatapoints;
[14806]85      [Storable]
[14788]86      public double finalMomentum;
[14806]87      [Storable]
[14788]88      public int momSwitchIter;
[14806]89      [Storable]
[14788]90      public int stopLyingIter;
[14806]91      [Storable]
[14788]92      public double theta;
[14806]93      [Storable]
[14788]94      public double eta;
[14806]95      [Storable]
[14788]96      public int newDimensions;
[14414]97
[14788]98      // for approximate version: sparse representation of similarity/distance matrix
[14806]99      [Storable]
[14788]100      public double[] valP; // similarity/distance
[14806]101      [Storable]
[14788]102      public int[] rowP; // row index
[14806]103      [Storable]
[14788]104      public int[] colP; // col index
[14414]105
[14788]106      // for exact version: dense representation of distance/similarity matrix
[14806]107      [Storable]
[14788]108      public double[,] p;
[14512]109
[14788]110      // mapped data
[14806]111      [Storable]
[14788]112      public double[,] newData;
[14414]113
[14806]114      [Storable]
[14788]115      public int iter;
[14806]116      [Storable]
[14788]117      public double currentMomentum;
[14414]118
[14788]119      // helper variables (updated in each iteration)
[14806]120      [Storable]
[14788]121      public double[,] gains;
[14806]122      [Storable]
[14788]123      public double[,] uY;
[14806]124      [Storable]
[14788]125      public double[,] dY;
[15207]126      #endregion
[14512]127
[15207]128      #region Constructors & Cloning
[14788]129      private TSNEState(TSNEState original, Cloner cloner) : base(original, cloner) {
[15207]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);
[14806]144        }
[15207]145        if (original.rowP != null) {
146          rowP = new int[original.rowP.Length];
147          Array.Copy(original.rowP, rowP, rowP.Length);
[14806]148        }
[15207]149        if (original.colP != null) {
150          colP = new int[original.colP.Length];
151          Array.Copy(original.colP, colP, colP.Length);
[14806]152        }
[15207]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);
[14806]156        }
[15207]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);
[14788]167      }
[14806]168
[14788]169      public override IDeepCloneable Clone(Cloner cloner) {
170        return new TSNEState(this, cloner);
171      }
[14414]172
[14807]173      [StorableConstructor]
[16565]174      private TSNEState(StorableConstructorFlag _) { }
[15532]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) {
[14788]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;
[15207]185        currentMomentum = momentum;
[14788]186        this.finalMomentum = finalMomentum;
187        this.eta = eta;
[14414]188
[14788]189        // initialize
[15532]190        noDatapoints = data.Count;
[15207]191        if (noDatapoints - 1 < 3 * perplexity)
[14806]192          throw new ArgumentException("Perplexity too large for the number of data points!");
[14788]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];
[15207]199        for (var i = 0; i < noDatapoints; i++)
[15532]200        for (var j = 0; j < newDimensions; j++)
201          gains[i, j] = 1.0;
[14788]202
203        p = null;
204        rowP = null;
205        colP = null;
206        valP = null;
207
208        //Calculate Similarities
[14858]209        if (exact) p = CalculateExactSimilarites(data, distance, perplexity);
[14788]210        else CalculateApproximateSimilarities(data, distance, perplexity, out rowP, out colP, out valP);
211
[14837]212        // Lie about the P-values (factor is 4 in the MATLAB implementation)
[15207]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;
[14788]215
216        // Initialize solution (randomly)
217        var rand = new NormalDistributedRandom(random, 0, 1);
[15207]218        for (var i = 0; i < noDatapoints; i++)
[15532]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        }
[14414]228      }
[15207]229      #endregion
[14414]230
[14788]231      public double EvaluateError() {
[15532]232        return exact ? EvaluateErrorExact(p, newData, noDatapoints, newDimensions) : EvaluateErrorApproximate(rowP, colP, valP, newData, theta);
[14788]233      }
[14512]234
[15207]235      #region Helpers
[15532]236      private static void CalculateApproximateSimilarities(IReadOnlyList<T> data, IDistance<T> distance, double perplexity, out int[] rowP, out int[] colP, out double[] valP) {
[14788]237        // Compute asymmetric pairwise input similarities
[15532]238        ComputeGaussianPerplexity(data, distance, out rowP, out colP, out valP, perplexity, (int) (3 * perplexity));
[14788]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;
[15532]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;
[14788]249      }
[15532]250      private static double[,] CalculateExactSimilarites(IReadOnlyList<T> data, IDistance<T> distance, double perplexity) {
[14788]251        // Compute similarities
[15532]252        var p = new double[data.Count, data.Count];
[14788]253        ComputeGaussianPerplexity(data, distance, p, perplexity);
254        // Symmetrize input similarities
[15532]255        for (var n = 0; n < data.Count; n++) {
256          for (var m = n + 1; m < data.Count; m++) {
[14788]257            p[n, m] += p[m, n];
258            p[m, n] = p[n, m];
259          }
260        }
261        var sumP = .0;
[15532]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        }
[14788]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) {
[15207]275        if (perplexity > k) throw new ArgumentException("Perplexity should be lower than k!");
[14512]276
[15207]277        var n = x.Count;
[14788]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;
[15207]284        for (var i = 0; i < n; i++) rowP[i + 1] = rowP[i] + k;
[14512]285
[14788]286        var objX = new List<IndexedItem<T>>();
[15207]287        for (var i = 0; i < n; i++) objX.Add(new IndexedItem<T>(i, x[i]));
[14512]288
[14788]289        // Build ball tree on data set
[14837]290        var tree = new VantagePointTree<IndexedItem<T>>(new IndexedItemDistance<T>(distance), objX);
[14742]291
[14788]292        // Loop over all points to find nearest neighbors
[15207]293        for (var i = 0; i < n; i++) {
[14788]294          IList<IndexedItem<T>> indices;
295          IList<double> distances;
[14742]296
[14788]297          // Find nearest neighbors
298          tree.Search(objX[i], k + 1, out indices, out distances);
[14512]299
[14788]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;
[15207]305          const double tol = 1e-5;
[14512]306
[14788]307          // Iterate until we found a good perplexity
[15532]308          var iter = 0;
309          double sumP = 0;
[15207]310          while (!found && iter < 200) {
[14788]311            // Compute Gaussian kernel row
[15207]312            for (var m = 0; m < k; m++) curP[m] = Math.Exp(-beta * distances[m + 1]);
[14512]313
[14788]314            // Compute entropy of current row
315            sumP = double.Epsilon;
[15207]316            for (var m = 0; m < k; m++) sumP += curP[m];
[14788]317            var h = .0;
[15207]318            for (var m = 0; m < k; m++) h += beta * (distances[m + 1] * curP[m]);
[14788]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);
[15207]323            if (hdiff < tol && -hdiff < tol) {
[14788]324              found = true;
[15532]325            }
326            else {
[15207]327              if (hdiff > 0) {
[14788]328                minBeta = beta;
[15207]329                if (maxBeta.IsAlmost(double.MaxValue) || maxBeta.IsAlmost(double.MinValue))
[14788]330                  beta *= 2.0;
331                else
332                  beta = (beta + maxBeta) / 2.0;
[15532]333              }
334              else {
[14788]335                maxBeta = beta;
[15207]336                if (minBeta.IsAlmost(double.MinValue) || minBeta.IsAlmost(double.MaxValue))
[14788]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
[15207]348          for (var m = 0; m < k; m++) curP[m] /= sumP;
349          for (var m = 0; m < k; m++) {
[14788]350            colP[rowP[i] + m] = indices[m + 1].Index;
351            valP[rowP[i] + m] = curP[m];
352          }
[14512]353        }
354      }
[15532]355      private static void ComputeGaussianPerplexity(IReadOnlyList<T> x, IDistance<T> distance, double[,] p, double perplexity) {
[14788]356        // Compute the distance matrix
357        var dd = ComputeDistances(x, distance);
358
[15532]359        var n = x.Count;
[14788]360        // Compute the Gaussian kernel row by row
[15207]361        for (var i = 0; i < n; i++) {
[14788]362          // Initialize some variables
363          var found = false;
364          var beta = 1.0;
[14837]365          var minBeta = double.MinValue;
[14788]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;
[15532]372          while (!found && iter < 200) { // 200 iterations as in tSNE implementation by van der Maarten
[14788]373
374            // Compute Gaussian kernel row
[15207]375            for (var m = 0; m < n; m++) p[i, m] = Math.Exp(-beta * dd[i][m]);
[14788]376            p[i, i] = double.Epsilon;
377
378            // Compute entropy of current row
379            sumP = double.Epsilon;
[15207]380            for (var m = 0; m < n; m++) sumP += p[i, m];
[14788]381            var h = 0.0;
[15207]382            for (var m = 0; m < n; m++) h += beta * (dd[i][m] * p[i, m]);
[14788]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);
[15207]387            if (hdiff < tol && -hdiff < tol) {
[14788]388              found = true;
[15532]389            }
390            else {
[15207]391              if (hdiff > 0) {
[14788]392                minBeta = beta;
[15207]393                if (maxBeta.IsAlmost(double.MaxValue) || maxBeta.IsAlmost(double.MinValue))
[14788]394                  beta *= 2.0;
395                else
396                  beta = (beta + maxBeta) / 2.0;
[15532]397              }
398              else {
[14788]399                maxBeta = beta;
[15207]400                if (minBeta.IsAlmost(double.MinValue) || minBeta.IsAlmost(double.MaxValue))
[14788]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
[15207]412          for (var m = 0; m < n; m++) p[i, m] /= sumP;
[14512]413        }
414      }
[15532]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];
[14806]419          // all distances must be symmetric
[15207]420          for (var c = 0; c < r; c++) {
[14806]421            rowV[c] = res[c][r];
422          }
423          rowV[r] = 0.0; // distance to self is zero for all distances
[15532]424          for (var c = r + 1; c < x.Count; c++) {
[14806]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();
[14788]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];
[14837]436        ComputeSquaredEuclideanDistance(y, n, d, dd);
[14414]437
[14788]438        // Compute Q-matrix and normalization sum
439        var sumQ = double.Epsilon;
[15207]440        for (var n1 = 0; n1 < n; n1++) {
441          for (var m = 0; m < n; m++) {
442            if (n1 != m) {
[14788]443              q[n1, m] = 1 / (1 + dd[n1, m]);
444              sumQ += q[n1, m];
[15532]445            }
446            else q[n1, m] = double.Epsilon;
[14788]447          }
448        }
[15207]449        for (var i = 0; i < n; i++) for (var j = 0; j < n; j++) q[i, j] /= sumQ;
[14414]450
[14788]451        // Sum t-SNE error
452        var c = .0;
[15207]453        for (var i = 0; i < n; i++)
[15532]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        }
[14788]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];
[15207]465        var sumQ = 0.0;
466        for (var i = 0; i < n; i++) tree.ComputeNonEdgeForces(i, theta, buff, ref sumQ);
[14414]467
[14788]468        // Loop over all edges to compute t-SNE error
469        var c = .0;
[15207]470        for (var k = 0; k < n; k++) {
471          for (var i = rowP[k]; i < rowP[k + 1]; i++) {
[14788]472            var q = .0;
[15207]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];
[14837]476            q = (1.0 / (1.0 + q)) / sumQ;
[14788]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];
[15207]486        for (var j = 0; j < n; j++) {
487          for (var i = rowP[j]; i < rowP[j + 1]; i++) {
[14788]488            // Check whether element (col_P[i], n) is present
489            var present = false;
[15207]490            for (var m = rowP[colP[i]]; m < rowP[colP[i] + 1]; m++) {
491              if (colP[m] == j) present = true;
[14788]492            }
[15207]493            if (present) rowCounts[j]++;
[14788]494            else {
495              rowCounts[j]++;
496              rowCounts[colP[i]]++;
497            }
498          }
499        }
500        var noElem = 0;
[15207]501        for (var i = 0; i < n; i++) noElem += rowCounts[i];
[14414]502
[14788]503        // Allocate memory for symmetrized matrix
504        symRowP = new int[n + 1];
505        symColP = new int[noElem];
506        symValP = new double[noElem];
[14414]507
[14788]508        // Construct new row indices for symmetric matrix
509        symRowP[0] = 0;
[15207]510        for (var i = 0; i < n; i++) symRowP[i + 1] = symRowP[i] + rowCounts[i];
[14788]511
512        // Fill the result matrix
513        var offset = new int[n];
[15207]514        for (var j = 0; j < n; j++) {
[15532]515          for (var i = rowP[j]; i < rowP[j + 1]; i++) { // considering element(n, colP[i])
[14788]516
517            // Check whether element (col_P[i], n) is present
518            var present = false;
[15207]519            for (var m = rowP[colP[i]]; m < rowP[colP[i] + 1]; m++) {
520              if (colP[m] != j) continue;
[14788]521              present = true;
[15207]522              if (j > colP[i]) continue; // make sure we do not add elements twice
[14788]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];
[14414]527            }
[14788]528
529            // If (colP[i], n) is not present, there is no addition involved
[15207]530            if (!present) {
[14788]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
[15207]538            if (present && (j > colP[i])) continue;
[14788]539            offset[j]++;
[15207]540            if (colP[i] != j) offset[colP[i]]++;
[14414]541          }
542        }
543
[15207]544        for (var i = 0; i < noElem; i++) symValP[i] /= 2.0;
[14414]545      }
[15207]546      #endregion
[14807]547    }
[14788]548
[14807]549    /// <summary>
[15207]550    /// Static interface to tSNE
[14807]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,
[15532]567      double theta = 0, int stopLyingIter = 0, int momSwitchIter = 0, double momentum = .5,
[15207]568      double finalMomentum = .8, double eta = 10.0
[15532]569    ) {
[14807]570      var state = CreateState(data, distance, random, newDimensions, perplexity,
571        theta, stopLyingIter, momSwitchIter, momentum, finalMomentum, eta);
572
[15207]573      for (var i = 0; i < iterations - 1; i++) {
[14807]574        Iterate(state);
575      }
576      return Iterate(state);
[14414]577    }
[14785]578
[14807]579    public static TSNEState CreateState(T[] data, IDistance<T> distance, IRandom random,
580      int newDimensions = 2, double perplexity = 25, double theta = 0,
[15207]581      int stopLyingIter = 0, int momSwitchIter = 0, double momentum = .5,
[15532]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);
[14788]585    }
[14414]586
[14788]587    public static double[,] Iterate(TSNEState state) {
[15207]588      if (state.exact)
[14788]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);
[14414]592
[14788]593      // Update gains
[15207]594      for (var i = 0; i < state.noDatapoints; i++) {
595        for (var j = 0; j < state.newDimensions; j++) {
[14788]596          state.gains[i, j] = Math.Sign(state.dY[i, j]) != Math.Sign(state.uY[i, j])
[15532]597            ? state.gains[i, j] + .2 // +0.2 nd *0.8 are used in two separate implementations of tSNE -> seems to be correct
[14837]598            : state.gains[i, j] * .8;
[15207]599          if (state.gains[i, j] < .01) state.gains[i, j] = .01;
[14414]600        }
[14788]601      }
[14414]602
[14788]603      // Perform gradient update (with momentum and gains)
[15207]604      for (var i = 0; i < state.noDatapoints; i++)
[15532]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];
[14788]607
[15207]608      for (var i = 0; i < state.noDatapoints; i++)
[15532]609      for (var j = 0; j < state.newDimensions; j++)
610        state.newData[i, j] = state.newData[i, j] + state.uY[i, j];
[14788]611
612      // Make solution zero-mean
613      ZeroMean(state.newData);
[14807]614
[14788]615      // Stop lying about the P-values after a while, and switch momentum
[15207]616      if (state.iter == state.stopLyingIter) {
617        if (state.exact)
618          for (var i = 0; i < state.noDatapoints; i++)
[15532]619          for (var j = 0; j < state.noDatapoints; j++)
620            state.p[i, j] /= 12.0;
[14788]621        else
[15207]622          for (var i = 0; i < state.rowP[state.noDatapoints]; i++)
[14837]623            state.valP[i] /= 12.0;
[14414]624      }
[14788]625
[15207]626      if (state.iter == state.momSwitchIter)
[14788]627        state.currentMomentum = state.finalMomentum;
628
629      state.iter++;
630      return state.newData;
[14414]631    }
[14785]632
[15207]633    #region Helpers
[14788]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);
[15207]636      var sumQ = 0.0;
[14788]637      var posF = new double[n, d];
638      var negF = new double[n, d];
[15207]639      SpacePartitioningTree.ComputeEdgeForces(rowP, colP, valP, n, posF, y, d);
[14788]640      var row = new double[d];
[15207]641      for (var n1 = 0; n1 < n; n1++) {
642        Array.Clear(row, 0, row.Length);
[14788]643        tree.ComputeNonEdgeForces(n1, theta, row, ref sumQ);
[15207]644        Buffer.BlockCopy(row, 0, negF, (sizeof(double) * n1 * d), d * sizeof(double));
[14788]645      }
646
647      // Compute final t-SNE gradient
[14856]648      for (var i = 0; i < n; i++)
[15532]649      for (var j = 0; j < d; j++) {
650        dC[i, j] = posF[i, j] - negF[i, j] / sumQ;
651      }
[14414]652    }
[14785]653
[14414]654    private static void ComputeExactGradient(double[,] p, double[,] y, int n, int d, double[,] dC) {
655      // Make sure the current gradient contains zeros
[15207]656      for (var i = 0; i < n; i++) for (var j = 0; j < d; j++) dC[i, j] = 0.0;
[14414]657
658      // Compute the squared Euclidean distance matrix
659      var dd = new double[n, n];
[14837]660      ComputeSquaredEuclideanDistance(y, n, d, dd);
[14414]661
662      // Compute Q-matrix and normalization sum
663      var q = new double[n, n];
664      var sumQ = .0;
[15207]665      for (var n1 = 0; n1 < n; n1++) {
666        for (var m = 0; m < n; m++) {
667          if (n1 == m) continue;
[14414]668          q[n1, m] = 1 / (1 + dd[n1, m]);
669          sumQ += q[n1, m];
670        }
671      }
672
673      // Perform the computation of the gradient
[15207]674      for (var n1 = 0; n1 < n; n1++) {
675        for (var m = 0; m < n; m++) {
676          if (n1 == m) continue;
[14414]677          var mult = (p[n1, m] - q[n1, m] / sumQ) * q[n1, m];
[15207]678          for (var d1 = 0; d1 < d; d1++) {
[14414]679            dC[n1, d1] += (y[n1, d1] - y[m, d1]) * mult;
680          }
681        }
682      }
683    }
[14788]684
[14414]685    private static void ComputeSquaredEuclideanDistance(double[,] x, int n, int d, double[,] dd) {
686      var dataSums = new double[n];
[15207]687      for (var i = 0; i < n; i++) {
688        for (var j = 0; j < d; j++) {
[14414]689          dataSums[i] += x[i, j] * x[i, j];
690        }
691      }
[15207]692      for (var i = 0; i < n; i++) {
693        for (var m = 0; m < n; m++) {
[14414]694          dd[i, m] = dataSums[i] + dataSums[m];
695        }
696      }
[15207]697      for (var i = 0; i < n; i++) {
[14414]698        dd[i, i] = 0.0;
[15207]699        for (var m = i + 1; m < n; m++) {
[14414]700          dd[i, m] = 0.0;
[15207]701          for (var j = 0; j < d; j++) {
[14414]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];
[15207]714      for (var i = 0; i < n; i++) {
715        for (var j = 0; j < d; j++) {
[14414]716          mean[j] += x[i, j];
717        }
718      }
[15207]719      for (var i = 0; i < d; i++) {
[14414]720        mean[i] /= n;
721      }
722      // Subtract data mean
[15207]723      for (var i = 0; i < n; i++) {
724        for (var j = 0; j < d; j++) {
[14414]725          x[i, j] -= mean[j];
726        }
727      }
728    }
[15207]729    #endregion
[14414]730  }
[15532]731}
Note: See TracBrowser for help on using the repository browser.