Free cookie consent management tool by TermsFeed Policy Generator

Ignore:
Timestamp:
07/15/17 10:29:40 (7 years ago)
Author:
gkronber
Message:

#2699,#2700
merged r14862, r14863, r14911, r14936, r15156, r15157, r15158, r15164, r15169, r15207:15209, r15225, r15227, r15234, r15248 from trunk to stable

Location:
stable
Files:
3 edited
1 copied

Legend:

Unmodified
Added
Removed

  • stable

  • stable/HeuristicLab.Algorithms.DataAnalysis

  • stable/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/TSNEStatic.cs

    r14862 r15249  
    5656using System;
    5757using System.Collections.Generic;
    58 using System.Linq;
    5958using HeuristicLab.Collections;
    6059using HeuristicLab.Common;
    6160using HeuristicLab.Core;
    62 using HeuristicLab.Optimization;
    6361using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
    6462using HeuristicLab.Random;
    … …  
    7068    [StorableClass]
    7169    public sealed class TSNEState : DeepCloneable {
     70      #region Storables
    7271      // initialized once
    7372      [Storable]
    … …  
    122121      [Storable]
    123122      public double[,] dY;
    124 
     123      #endregion
     124
     125      #region Constructors & Cloning
    125126      private TSNEState(TSNEState original, Cloner cloner) : base(original, cloner) {
    126         this.distance = cloner.Clone(original.distance);
    127         this.random = cloner.Clone(original.random);
    128         this.perplexity = original.perplexity;
    129         this.exact = original.exact;
    130         this.noDatapoints = original.noDatapoints;
    131         this.finalMomentum = original.finalMomentum;
    132         this.momSwitchIter = original.momSwitchIter;
    133         this.stopLyingIter = original.stopLyingIter;
    134         this.theta = original.theta;
    135         this.eta = original.eta;
    136         this.newDimensions = original.newDimensions;
    137         if(original.valP != null) {
    138           this.valP = new double[original.valP.Length];
    139           Array.Copy(original.valP, this.valP, this.valP.Length);
    140         }
    141         if(original.rowP != null) {
    142           this.rowP = new int[original.rowP.Length];
    143           Array.Copy(original.rowP, this.rowP, this.rowP.Length);
    144         }
    145         if(original.colP != null) {
    146           this.colP = new int[original.colP.Length];
    147           Array.Copy(original.colP, this.colP, this.colP.Length);
    148         }
    149         if(original.p != null) {
    150           this.p = new double[original.p.GetLength(0), original.p.GetLength(1)];
    151           Array.Copy(original.p, this.p, this.p.Length);
    152         }
    153         this.newData = new double[original.newData.GetLength(0), original.newData.GetLength(1)];
    154         Array.Copy(original.newData, this.newData, this.newData.Length);
    155         this.iter = original.iter;
    156         this.currentMomentum = original.currentMomentum;
    157         this.gains = new double[original.gains.GetLength(0), original.gains.GetLength(1)];
    158         Array.Copy(original.gains, this.gains, this.gains.Length);
    159         this.uY = new double[original.uY.GetLength(0), original.uY.GetLength(1)];
    160         Array.Copy(original.uY, this.uY, this.uY.Length);
    161         this.dY = new double[original.dY.GetLength(0), original.dY.GetLength(1)];
    162         Array.Copy(original.dY, this.dY, this.dY.Length);
     127        distance = cloner.Clone(original.distance);
     128        random = cloner.Clone(original.random);
     129        perplexity = original.perplexity;
     130        exact = original.exact;
     131        noDatapoints = original.noDatapoints;
     132        finalMomentum = original.finalMomentum;
     133        momSwitchIter = original.momSwitchIter;
     134        stopLyingIter = original.stopLyingIter;
     135        theta = original.theta;
     136        eta = original.eta;
     137        newDimensions = original.newDimensions;
     138        if (original.valP != null) {
     139          valP = new double[original.valP.Length];
     140          Array.Copy(original.valP, valP, valP.Length);
     141        }
     142        if (original.rowP != null) {
     143          rowP = new int[original.rowP.Length];
     144          Array.Copy(original.rowP, rowP, rowP.Length);
     145        }
     146        if (original.colP != null) {
     147          colP = new int[original.colP.Length];
     148          Array.Copy(original.colP, colP, colP.Length);
     149        }
     150        if (original.p != null) {
     151          p = new double[original.p.GetLength(0), original.p.GetLength(1)];
     152          Array.Copy(original.p, p, p.Length);
     153        }
     154        newData = new double[original.newData.GetLength(0), original.newData.GetLength(1)];
     155        Array.Copy(original.newData, newData, newData.Length);
     156        iter = original.iter;
     157        currentMomentum = original.currentMomentum;
     158        gains = new double[original.gains.GetLength(0), original.gains.GetLength(1)];
     159        Array.Copy(original.gains, gains, gains.Length);
     160        uY = new double[original.uY.GetLength(0), original.uY.GetLength(1)];
     161        Array.Copy(original.uY, uY, uY.Length);
     162        dY = new double[original.dY.GetLength(0), original.dY.GetLength(1)];
     163        Array.Copy(original.dY, dY, dY.Length);
    163164      }
    164165
    … …  
    177178        this.stopLyingIter = stopLyingIter;
    178179        this.momSwitchIter = momSwitchIter;
    179         this.currentMomentum = momentum;
     180        currentMomentum = momentum;
    180181        this.finalMomentum = finalMomentum;
    181182        this.eta = eta;
    182183
    183 
    184184        // initialize
    185185        noDatapoints = data.Length;
    186         if(noDatapoints - 1 < 3 * perplexity)
     186        if (noDatapoints - 1 < 3 * perplexity)
    187187          throw new ArgumentException("Perplexity too large for the number of data points!");
    188188
    … …  
    192192        uY = new double[noDatapoints, newDimensions];
    193193        gains = new double[noDatapoints, newDimensions];
    194         for(var i = 0; i < noDatapoints; i++)
    195           for(var j = 0; j < newDimensions; j++)
     194        for (var i = 0; i < noDatapoints; i++)
     195          for (var j = 0; j < newDimensions; j++)
    196196            gains[i, j] = 1.0;
    197197
    … …  
    206206
    207207        // Lie about the P-values (factor is 4 in the MATLAB implementation)
    208         if(exact) for(var i = 0; i < noDatapoints; i++) for(var j = 0; j < noDatapoints; j++) p[i, j] *= 12.0;
    209         else for(var i = 0; i < rowP[noDatapoints]; i++) valP[i] *= 12.0;
     208        if (exact) for (var i = 0; i < noDatapoints; i++) for (var j = 0; j < noDatapoints; j++) p[i, j] *= 12.0;
     209        else for (var i = 0; i < rowP[noDatapoints]; i++) valP[i] *= 12.0;
    210210
    211211        // Initialize solution (randomly)
    212212        var rand = new NormalDistributedRandom(random, 0, 1);
    213         for(var i = 0; i < noDatapoints; i++)
    214           for(var j = 0; j < newDimensions; j++)
     213        for (var i = 0; i < noDatapoints; i++)
     214          for (var j = 0; j < newDimensions; j++)
    215215            newData[i, j] = rand.NextDouble() * .0001;
    216216      }
     217      #endregion
    217218
    218219      public double EvaluateError() {
    … …  
    222223      }
    223224
     225      #region Helpers
    224226      private static void CalculateApproximateSimilarities(T[] data, IDistance<T> distance, double perplexity, out int[] rowP, out int[] colP, out double[] valP) {
    225227        // Compute asymmetric pairwise input similarities
    … …  
    233235        valP = sValP;
    234236        var sumP = .0;
    235         for(var i = 0; i < rowP[data.Length]; i++) sumP += valP[i];
    236         for(var i = 0; i < rowP[data.Length]; i++) valP[i] /= sumP;
     237        for (var i = 0; i < rowP[data.Length]; i++) sumP += valP[i];
     238        for (var i = 0; i < rowP[data.Length]; i++) valP[i] /= sumP;
    237239      }
    238240
    … …  
    242244        ComputeGaussianPerplexity(data, distance, p, perplexity);
    243245        // Symmetrize input similarities
    244         for(var n = 0; n < data.Length; n++) {
    245           for(var m = n + 1; m < data.Length; m++) {
     246        for (var n = 0; n < data.Length; n++) {
     247          for (var m = n + 1; m < data.Length; m++) {
    246248            p[n, m] += p[m, n];
    247249            p[m, n] = p[n, m];
    … …  
    249251        }
    250252        var sumP = .0;
    251         for(var i = 0; i < data.Length; i++) for(var j = 0; j < data.Length; j++) sumP += p[i, j];
    252         for(var i = 0; i < data.Length; i++) for(var j = 0; j < data.Length; j++) p[i, j] /= sumP;
     253        for (var i = 0; i < data.Length; i++) for (var j = 0; j < data.Length; j++) sumP += p[i, j];
     254        for (var i = 0; i < data.Length; i++) for (var j = 0; j < data.Length; j++) p[i, j] /= sumP;
    253255        return p;
    254256      }
    255257
    256258      private static void ComputeGaussianPerplexity(IReadOnlyList<T> x, IDistance<T> distance, out int[] rowP, out int[] colP, out double[] valP, double perplexity, int k) {
    257         if(perplexity > k) throw new ArgumentException("Perplexity should be lower than k!");
    258 
    259         int n = x.Count;
     259        if (perplexity > k) throw new ArgumentException("Perplexity should be lower than k!");
     260
     261        var n = x.Count;
    260262        // Allocate the memory we need
    261263        rowP = new int[n + 1];
    … …  
    264266        var curP = new double[n - 1];
    265267        rowP[0] = 0;
    266         for(var i = 0; i < n; i++) rowP[i + 1] = rowP[i] + k;
     268        for (var i = 0; i < n; i++) rowP[i + 1] = rowP[i] + k;
    267269
    268270        var objX = new List<IndexedItem<T>>();
    269         for(var i = 0; i < n; i++) objX.Add(new IndexedItem<T>(i, x[i]));
     271        for (var i = 0; i < n; i++) objX.Add(new IndexedItem<T>(i, x[i]));
    270272
    271273        // Build ball tree on data set
    … …  
    273275
    274276        // Loop over all points to find nearest neighbors
    275         for(var i = 0; i < n; i++) {
     277        for (var i = 0; i < n; i++) {
    276278          IList<IndexedItem<T>> indices;
    277279          IList<double> distances;
    … …  
    285287          var minBeta = double.MinValue;
    286288          var maxBeta = double.MaxValue;
    287           const double tol = 1e-5; 
     289          const double tol = 1e-5;
    288290
    289291          // Iterate until we found a good perplexity
    290292          var iter = 0; double sumP = 0;
    291           while(!found && iter < 200) { 
     293          while (!found && iter < 200) {
    292294
    293295            // Compute Gaussian kernel row
    294             for(var m = 0; m < k; m++) curP[m] = Math.Exp(-beta * distances[m + 1]);
     296            for (var m = 0; m < k; m++) curP[m] = Math.Exp(-beta * distances[m + 1]);
    295297
    296298            // Compute entropy of current row
    297299            sumP = double.Epsilon;
    298             for(var m = 0; m < k; m++) sumP += curP[m];
     300            for (var m = 0; m < k; m++) sumP += curP[m];
    299301            var h = .0;
    300             for(var m = 0; m < k; m++) h += beta * (distances[m + 1] * curP[m]);
     302            for (var m = 0; m < k; m++) h += beta * (distances[m + 1] * curP[m]);
    301303            h = h / sumP + Math.Log(sumP);
    302304
    303305            // Evaluate whether the entropy is within the tolerance level
    304306            var hdiff = h - Math.Log(perplexity);
    305             if(hdiff < tol && -hdiff < tol) {
     307            if (hdiff < tol && -hdiff < tol) {
    306308              found = true;
    307309            } else {
    308               if(hdiff > 0) {
     310              if (hdiff > 0) {
    309311                minBeta = beta;
    310                 if(maxBeta.IsAlmost(double.MaxValue) || maxBeta.IsAlmost(double.MinValue))
     312                if (maxBeta.IsAlmost(double.MaxValue) || maxBeta.IsAlmost(double.MinValue))
    311313                  beta *= 2.0;
    312314                else
    … …  
    314316              } else {
    315317                maxBeta = beta;
    316                 if(minBeta.IsAlmost(double.MinValue) || minBeta.IsAlmost(double.MaxValue))
     318                if (minBeta.IsAlmost(double.MinValue) || minBeta.IsAlmost(double.MaxValue))
    317319                  beta /= 2.0;
    318320                else
    … …  
    326328
    327329          // Row-normalize current row of P and store in matrix
    328           for(var m = 0; m < k; m++) curP[m] /= sumP;
    329           for(var m = 0; m < k; m++) {
     330          for (var m = 0; m < k; m++) curP[m] /= sumP;
     331          for (var m = 0; m < k; m++) {
    330332            colP[rowP[i] + m] = indices[m + 1].Index;
    331333            valP[rowP[i] + m] = curP[m];
    … …  
    337339        var dd = ComputeDistances(x, distance);
    338340
    339         int n = x.Length;
     341        var n = x.Length;
    340342        // Compute the Gaussian kernel row by row
    341         for(var i = 0; i < n; i++) {
     343        for (var i = 0; i < n; i++) {
    342344          // Initialize some variables
    343345          var found = false;
    … …  
    350352          // Iterate until we found a good perplexity
    351353          var iter = 0;
    352           while(!found && iter < 200) {      // 200 iterations as in tSNE implementation by van der Maarten
     354          while (!found && iter < 200) {      // 200 iterations as in tSNE implementation by van der Maarten
    353355
    354356            // Compute Gaussian kernel row
    355             for(var m = 0; m < n; m++) p[i, m] = Math.Exp(-beta * dd[i][m]);
     357            for (var m = 0; m < n; m++) p[i, m] = Math.Exp(-beta * dd[i][m]);
    356358            p[i, i] = double.Epsilon;
    357359
    358360            // Compute entropy of current row
    359361            sumP = double.Epsilon;
    360             for(var m = 0; m < n; m++) sumP += p[i, m];
     362            for (var m = 0; m < n; m++) sumP += p[i, m];
    361363            var h = 0.0;
    362             for(var m = 0; m < n; m++) h += beta * (dd[i][m] * p[i, m]);
     364            for (var m = 0; m < n; m++) h += beta * (dd[i][m] * p[i, m]);
    363365            h = h / sumP + Math.Log(sumP);
    364366
    365367            // Evaluate whether the entropy is within the tolerance level
    366368            var hdiff = h - Math.Log(perplexity);
    367             if(hdiff < tol && -hdiff < tol) {
     369            if (hdiff < tol && -hdiff < tol) {
    368370              found = true;
    369371            } else {
    370               if(hdiff > 0) {
     372              if (hdiff > 0) {
    371373                minBeta = beta;
    372                 if(maxBeta.IsAlmost(double.MaxValue) || maxBeta.IsAlmost(double.MinValue))
     374                if (maxBeta.IsAlmost(double.MaxValue) || maxBeta.IsAlmost(double.MinValue))
    373375                  beta *= 2.0;
    374376                else
    … …  
    376378              } else {
    377379                maxBeta = beta;
    378                 if(minBeta.IsAlmost(double.MinValue) || minBeta.IsAlmost(double.MaxValue))
     380                if (minBeta.IsAlmost(double.MinValue) || minBeta.IsAlmost(double.MaxValue))
    379381                  beta /= 2.0;
    380382                else
    … …  
    388390
    389391          // Row normalize P
    390           for(var m = 0; m < n; m++) p[i, m] /= sumP;
     392          for (var m = 0; m < n; m++) p[i, m] /= sumP;
    391393        }
    392394      }
    … …  
    394396      private static double[][] ComputeDistances(T[] x, IDistance<T> distance) {
    395397        var res = new double[x.Length][];
    396         for(int r = 0; r < x.Length; r++) {
     398        for (var r = 0; r < x.Length; r++) {
    397399          var rowV = new double[x.Length];
    398400          // all distances must be symmetric
    399           for(int c = 0; c < r; c++) {
     401          for (var c = 0; c < r; c++) {
    400402            rowV[c] = res[c][r];
    401403          }
    402404          rowV[r] = 0.0; // distance to self is zero for all distances
    403           for(int c = r + 1; c < x.Length; c++) {
     405          for (var c = r + 1; c < x.Length; c++) {
    404406            rowV[c] = distance.Get(x[r], x[c]);
    405407          }
    … …  
    418420        // Compute Q-matrix and normalization sum
    419421        var sumQ = double.Epsilon;
    420         for(var n1 = 0; n1 < n; n1++) {
    421           for(var m = 0; m < n; m++) {
    422             if(n1 != m) {
     422        for (var n1 = 0; n1 < n; n1++) {
     423          for (var m = 0; m < n; m++) {
     424            if (n1 != m) {
    423425              q[n1, m] = 1 / (1 + dd[n1, m]);
    424426              sumQ += q[n1, m];
    … …  
    426428          }
    427429        }
    428         for(var i = 0; i < n; i++) for(var j = 0; j < n; j++) q[i, j] /= sumQ;
     430        for (var i = 0; i < n; i++) for (var j = 0; j < n; j++) q[i, j] /= sumQ;
    429431
    430432        // Sum t-SNE error
    431433        var c = .0;
    432         for(var i = 0; i < n; i++)
    433           for(var j = 0; j < n; j++) {
     434        for (var i = 0; i < n; i++)
     435          for (var j = 0; j < n; j++) {
    434436            c += p[i, j] * Math.Log((p[i, j] + float.Epsilon) / (q[i, j] + float.Epsilon));
    435437          }
    … …  
    437439      }
    438440
    439       // The mapping of the approximate tSNE looks good but the error curve never changes.
    440441      private static double EvaluateErrorApproximate(IReadOnlyList<int> rowP, IReadOnlyList<int> colP, IReadOnlyList<double> valP, double[,] y, double theta) {
    441442        // Get estimate of normalization term
    … …  
    444445        var tree = new SpacePartitioningTree(y);
    445446        var buff = new double[d];
    446         double sumQ = 0.0;
    447         for(var i = 0; i < n; i++) tree.ComputeNonEdgeForces(i, theta, buff, ref sumQ);
     447        var sumQ = 0.0;
     448        for (var i = 0; i < n; i++) tree.ComputeNonEdgeForces(i, theta, buff, ref sumQ);
    448449
    449450        // Loop over all edges to compute t-SNE error
    450451        var c = .0;
    451         for(var k = 0; k < n; k++) {
    452           for(var i = rowP[k]; i < rowP[k + 1]; i++) {
     452        for (var k = 0; k < n; k++) {
     453          for (var i = rowP[k]; i < rowP[k + 1]; i++) {
    453454            var q = .0;
    454             for(var j = 0; j < d; j++) buff[j] = y[k, j];
    455             for(var j = 0; j < d; j++) buff[j] -= y[colP[i], j];
    456             for(var j = 0; j < d; j++) q += buff[j] * buff[j];
     455            for (var j = 0; j < d; j++) buff[j] = y[k, j];
     456            for (var j = 0; j < d; j++) buff[j] -= y[colP[i], j];
     457            for (var j = 0; j < d; j++) q += buff[j] * buff[j];
    457458            q = (1.0 / (1.0 + q)) / sumQ;
    458459            c += valP[i] * Math.Log((valP[i] + float.Epsilon) / (q + float.Epsilon));
    … …  
    466467        var n = rowP.Count - 1;
    467468        var rowCounts = new int[n];
    468         for(var j = 0; j < n; j++) {
    469           for(var i = rowP[j]; i < rowP[j + 1]; i++) {
     469        for (var j = 0; j < n; j++) {
     470          for (var i = rowP[j]; i < rowP[j + 1]; i++) {
    470471
    471472            // Check whether element (col_P[i], n) is present
    472473            var present = false;
    473             for(var m = rowP[colP[i]]; m < rowP[colP[i] + 1]; m++) {
    474               if(colP[m] == j) present = true;
     474            for (var m = rowP[colP[i]]; m < rowP[colP[i] + 1]; m++) {
     475              if (colP[m] == j) present = true;
    475476            }
    476             if(present) rowCounts[j]++;
     477            if (present) rowCounts[j]++;
    477478            else {
    478479              rowCounts[j]++;
    … …  
    482483        }
    483484        var noElem = 0;
    484         for(var i = 0; i < n; i++) noElem += rowCounts[i];
     485        for (var i = 0; i < n; i++) noElem += rowCounts[i];
    485486
    486487        // Allocate memory for symmetrized matrix
    … …  
    491492        // Construct new row indices for symmetric matrix
    492493        symRowP[0] = 0;
    493         for(var i = 0; i < n; i++) symRowP[i + 1] = symRowP[i] + rowCounts[i];
     494        for (var i = 0; i < n; i++) symRowP[i + 1] = symRowP[i] + rowCounts[i];
    494495
    495496        // Fill the result matrix
    496497        var offset = new int[n];
    497         for(var j = 0; j < n; j++) {
    498           for(var i = rowP[j]; i < rowP[j + 1]; i++) {                                  // considering element(n, colP[i])
     498        for (var j = 0; j < n; j++) {
     499          for (var i = rowP[j]; i < rowP[j + 1]; i++) {                                  // considering element(n, colP[i])
    499500
    500501            // Check whether element (col_P[i], n) is present
    501502            var present = false;
    502             for(var m = rowP[colP[i]]; m < rowP[colP[i] + 1]; m++) {
    503               if(colP[m] != j) continue;
     503            for (var m = rowP[colP[i]]; m < rowP[colP[i] + 1]; m++) {
     504              if (colP[m] != j) continue;
    504505              present = true;
    505               if(j > colP[i]) continue; // make sure we do not add elements twice
     506              if (j > colP[i]) continue; // make sure we do not add elements twice
    506507              symColP[symRowP[j] + offset[j]] = colP[i];
    507508              symColP[symRowP[colP[i]] + offset[colP[i]]] = j;
    … …  
    511512
    512513            // If (colP[i], n) is not present, there is no addition involved
    513             if(!present) {
     514            if (!present) {
    514515              symColP[symRowP[j] + offset[j]] = colP[i];
    515516              symColP[symRowP[colP[i]] + offset[colP[i]]] = j;
    … …  
    519520
    520521            // Update offsets
    521             if(present && (j > colP[i])) continue;
     522            if (present && (j > colP[i])) continue;
    522523            offset[j]++;
    523             if(colP[i] != j) offset[colP[i]]++;
    524           }
    525         }
    526 
    527         for(var i = 0; i < noElem; i++) symValP[i] /= 2.0;
    528       }
     524            if (colP[i] != j) offset[colP[i]]++;
     525          }
     526        }
     527
     528        for (var i = 0; i < noElem; i++) symValP[i] /= 2.0;
     529      }
     530      #endregion
    529531    }
    530532
    531533    /// <summary>
    532     /// Simple interface to tSNE
     534    /// Static interface to tSNE
    533535    /// </summary>
    534536    /// <param name="data"></param>
    … …  
    548550      int newDimensions = 2, double perplexity = 25, int iterations = 1000,
    549551      double theta = 0,
    550       int stopLyingIter = 250, int momSwitchIter = 250, double momentum = .5,
    551       double finalMomentum = .8, double eta = 200.0
     552      int stopLyingIter = 0, int momSwitchIter = 0, double momentum = .5,
     553      double finalMomentum = .8, double eta = 10.0
    552554      ) {
    553555      var state = CreateState(data, distance, random, newDimensions, perplexity,
    554556        theta, stopLyingIter, momSwitchIter, momentum, finalMomentum, eta);
    555557
    556       for(int i = 0; i < iterations - 1; i++) {
     558      for (var i = 0; i < iterations - 1; i++) {
    557559        Iterate(state);
    558560      }
    … …  
    562564    public static TSNEState CreateState(T[] data, IDistance<T> distance, IRandom random,
    563565      int newDimensions = 2, double perplexity = 25, double theta = 0,
    564       int stopLyingIter = 250, int momSwitchIter = 250, double momentum = .5,
    565       double finalMomentum = .8, double eta = 200.0
     566      int stopLyingIter = 0, int momSwitchIter = 0, double momentum = .5,
     567      double finalMomentum = .8, double eta = 10.0
    566568      ) {
    567569      return new TSNEState(data, distance, random, newDimensions, perplexity, theta, stopLyingIter, momSwitchIter, momentum, finalMomentum, eta);
    568570    }
    569571
    570 
    571572    public static double[,] Iterate(TSNEState state) {
    572       if(state.exact)
     573      if (state.exact)
    573574        ComputeExactGradient(state.p, state.newData, state.noDatapoints, state.newDimensions, state.dY);
    574575      else
    … …  
    576577
    577578      // Update gains
    578       for(var i = 0; i < state.noDatapoints; i++) {
    579         for(var j = 0; j < state.newDimensions; j++) {
     579      for (var i = 0; i < state.noDatapoints; i++) {
     580        for (var j = 0; j < state.newDimensions; j++) {
    580581          state.gains[i, j] = Math.Sign(state.dY[i, j]) != Math.Sign(state.uY[i, j])
    581582            ? state.gains[i, j] + .2  // +0.2 nd *0.8 are used in two separate implementations of tSNE -> seems to be correct
    582583            : state.gains[i, j] * .8;
    583584
    584           if(state.gains[i, j] < .01) state.gains[i, j] = .01;
     585          if (state.gains[i, j] < .01) state.gains[i, j] = .01;
    585586        }
    586587      }
    … …  
    588589
    589590      // Perform gradient update (with momentum and gains)
    590       for(var i = 0; i < state.noDatapoints; i++)
    591         for(var j = 0; j < state.newDimensions; j++)
     591      for (var i = 0; i < state.noDatapoints; i++)
     592        for (var j = 0; j < state.newDimensions; j++)
    592593          state.uY[i, j] = state.currentMomentum * state.uY[i, j] - state.eta * state.gains[i, j] * state.dY[i, j];
    593594
    594       for(var i = 0; i < state.noDatapoints; i++)
    595         for(var j = 0; j < state.newDimensions; j++)
     595      for (var i = 0; i < state.noDatapoints; i++)
     596        for (var j = 0; j < state.newDimensions; j++)
    596597          state.newData[i, j] = state.newData[i, j] + state.uY[i, j];
    597598
    … …  
    600601
    601602      // Stop lying about the P-values after a while, and switch momentum
    602       if(state.iter == state.stopLyingIter) {
    603         if(state.exact)
    604           for(var i = 0; i < state.noDatapoints; i++)
    605             for(var j = 0; j < state.noDatapoints; j++)
     603      if (state.iter == state.stopLyingIter) {
     604        if (state.exact)
     605          for (var i = 0; i < state.noDatapoints; i++)
     606            for (var j = 0; j < state.noDatapoints; j++)
    606607              state.p[i, j] /= 12.0;
    607608        else
    608           for(var i = 0; i < state.rowP[state.noDatapoints]; i++)
     609          for (var i = 0; i < state.rowP[state.noDatapoints]; i++)
    609610            state.valP[i] /= 12.0;
    610611      }
    611612
    612       if(state.iter == state.momSwitchIter)
     613      if (state.iter == state.momSwitchIter)
    613614        state.currentMomentum = state.finalMomentum;
    614615
    … …  
    617618    }
    618619
    619 
    620 
     620    #region Helpers
    621621    private static void ComputeApproximateGradient(int[] rowP, int[] colP, double[] valP, double[,] y, int n, int d, double[,] dC, double theta) {
    622622      var tree = new SpacePartitioningTree(y);
    623       double sumQ = 0.0;
     623      var sumQ = 0.0;
    624624      var posF = new double[n, d];
    625625      var negF = new double[n, d];
    626       tree.ComputeEdgeForces(rowP, colP, valP, n, posF);
     626      SpacePartitioningTree.ComputeEdgeForces(rowP, colP, valP, n, posF, y, d);
    627627      var row = new double[d];
    628       for(var n1 = 0; n1 < n; n1++) {
    629         Array.Clear(row,0,row.Length);
     628      for (var n1 = 0; n1 < n; n1++) {
     629        Array.Clear(row, 0, row.Length);
    630630        tree.ComputeNonEdgeForces(n1, theta, row, ref sumQ);
    631         Buffer.BlockCopy(row, 0, negF, (sizeof(double) * n1 * d), d*sizeof(double));
     631        Buffer.BlockCopy(row, 0, negF, (sizeof(double) * n1 * d), d * sizeof(double));
    632632      }
    633633
    634634      // Compute final t-SNE gradient
    635635      for (var i = 0; i < n; i++)
    636         for(var j = 0; j < d; j++) {
     636        for (var j = 0; j < d; j++) {
    637637          dC[i, j] = posF[i, j] - negF[i, j] / sumQ;
    638638        }
    … …  
    640640
    641641    private static void ComputeExactGradient(double[,] p, double[,] y, int n, int d, double[,] dC) {
    642 
    643642      // Make sure the current gradient contains zeros
    644       for(var i = 0; i < n; i++) for(var j = 0; j < d; j++) dC[i, j] = 0.0;
     643      for (var i = 0; i < n; i++) for (var j = 0; j < d; j++) dC[i, j] = 0.0;
    645644
    646645      // Compute the squared Euclidean distance matrix
    … …  
    651650      var q = new double[n, n];
    652651      var sumQ = .0;
    653       for(var n1 = 0; n1 < n; n1++) {
    654         for(var m = 0; m < n; m++) {
    655           if(n1 == m) continue;
     652      for (var n1 = 0; n1 < n; n1++) {
     653        for (var m = 0; m < n; m++) {
     654          if (n1 == m) continue;
    656655          q[n1, m] = 1 / (1 + dd[n1, m]);
    657656          sumQ += q[n1, m];
    … …  
    660659
    661660      // Perform the computation of the gradient
    662       for(var n1 = 0; n1 < n; n1++) {
    663         for(var m = 0; m < n; m++) {
    664           if(n1 == m) continue;
     661      for (var n1 = 0; n1 < n; n1++) {
     662        for (var m = 0; m < n; m++) {
     663          if (n1 == m) continue;
    665664          var mult = (p[n1, m] - q[n1, m] / sumQ) * q[n1, m];
    666           for(var d1 = 0; d1 < d; d1++) {
     665          for (var d1 = 0; d1 < d; d1++) {
    667666            dC[n1, d1] += (y[n1, d1] - y[m, d1]) * mult;
    668667          }
    … …  
    673672    private static void ComputeSquaredEuclideanDistance(double[,] x, int n, int d, double[,] dd) {
    674673      var dataSums = new double[n];
    675       for(var i = 0; i < n; i++) {
    676         for(var j = 0; j < d; j++) {
     674      for (var i = 0; i < n; i++) {
     675        for (var j = 0; j < d; j++) {
    677676          dataSums[i] += x[i, j] * x[i, j];
    678677        }
    679678      }
    680       for(var i = 0; i < n; i++) {
    681         for(var m = 0; m < n; m++) {
     679      for (var i = 0; i < n; i++) {
     680        for (var m = 0; m < n; m++) {
    682681          dd[i, m] = dataSums[i] + dataSums[m];
    683682        }
    684683      }
    685       for(var i = 0; i < n; i++) {
     684      for (var i = 0; i < n; i++) {
    686685        dd[i, i] = 0.0;
    687         for(var m = i + 1; m < n; m++) {
     686        for (var m = i + 1; m < n; m++) {
    688687          dd[i, m] = 0.0;
    689           for(var j = 0; j < d; j++) {
     688          for (var j = 0; j < d; j++) {
    690689            dd[i, m] += (x[i, j] - x[m, j]) * (x[i, j] - x[m, j]);
    691690          }
    … …  
    700699      var d = x.GetLength(1);
    701700      var mean = new double[d];
    702       for(var i = 0; i < n; i++) {
    703         for(var j = 0; j < d; j++) {
     701      for (var i = 0; i < n; i++) {
     702        for (var j = 0; j < d; j++) {
    704703          mean[j] += x[i, j];
    705704        }
    706705      }
    707       for(var i = 0; i < d; i++) {
     706      for (var i = 0; i < d; i++) {
    708707        mean[i] /= n;
    709708      }
    710709      // Subtract data mean
    711       for(var i = 0; i < n; i++) {
    712         for(var j = 0; j < d; j++) {
     710      for (var i = 0; i < n; i++) {
     711        for (var j = 0; j < d; j++) {
    713712          x[i, j] -= mean[j];
    714713        }
    715714      }
    716715    }
     716    #endregion
    717717  }
    718718}
Note: See TracChangeset for help on using the changeset viewer.