Changeset 14806


Ignore:
Timestamp:
03/30/17 17:54:45 (5 years ago)
Author:
gkronber
Message:

#2700: worked on tSNE, storable and cloning for tSNE state. Added some TODO comments while reviewing.

Location:
branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/Interfaces/TSNEInterfaces/ISpacePartitioningTree.cs

    r14788 r14806  
    2525  public interface ISpacePartitioningTree {
    2626
    27     double[,] Data { set; }
    2827    int GetDepth();
    29     ISpacePartitioningTree GetParent();
    3028    bool Insert(int newIndex);
    3129    void Subdivide();
  • branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/SpacePartitioningTree.cs

    r14788 r14806  
    7777
    7878    // Indices in this space-partitioning tree node, corresponding center-of-mass, and list of all children
    79     public double[,] Data { private get; set; } // TODO public setter, private getter?
     79    private double[,] data;
    8080
    8181    private double[] centerOfMass;
     
    9191      var meanY = new double[d];
    9292      var minY = new double[d];
    93       for (var i = 0; i < d; i++) minY[i] = double.MaxValue;
     93      for(var i = 0; i < d; i++) minY[i] = double.MaxValue;
    9494      var maxY = new double[d];
    95       for (var i = 0; i < d; i++) maxY[i] = double.MinValue;
    96       for (uint i = 0; i < n; i++) {
    97         for (uint j = 0; j < d; j++) {
     95      for(var i = 0; i < d; i++) maxY[i] = double.MinValue;
     96      for(uint i = 0; i < n; i++) {
     97        for(uint j = 0; j < d; j++) {
    9898          meanY[j] = inpData[i, j];
    99           if (inpData[i, j] < minY[j]) minY[j] = inpData[i, j];
    100           if (inpData[i, j] > maxY[j]) maxY[j] = inpData[i, j];
    101         }
    102       }
    103       for (var i = 0; i < d; i++) meanY[i] /= n;
     99          if(inpData[i, j] < minY[j]) minY[j] = inpData[i, j];
     100          if(inpData[i, j] > maxY[j]) maxY[j] = inpData[i, j];
     101        }
     102      }
     103      for(var i = 0; i < d; i++) meanY[i] /= n;
    104104      var width = new double[d];
    105       for (var i = 0; i < d; i++) width[i] = Math.Max(maxY[i] - meanY[i], meanY[i] - maxY[i]) + 1e-5;
     105      for(var i = 0; i < d; i++) width[i] = Math.Max(maxY[i] - meanY[i], meanY[i] - maxY[i]) + 1e-5; // TODO constant?
    106106      Init(null, inpData, meanY, width);
    107107      Fill(n);
     
    122122      // Ignore objects which do not belong in this quad tree
    123123      var point = new double[dimension];
    124       Buffer.BlockCopy(Data, sizeof(double) * dimension * newIndex, point, 0, sizeof(double) * dimension);
    125       if (!boundary.ContainsPoint(point)) return false;
     124      Buffer.BlockCopy(data, sizeof(double) * dimension * newIndex, point, 0, sizeof(double) * dimension);
     125      if(!boundary.ContainsPoint(point)) return false;
    126126      cumulativeSize++;
    127127      // Online update of cumulative size and center-of-mass
    128128      var mult1 = (double)(cumulativeSize - 1) / cumulativeSize;
    129129      var mult2 = 1.0 / cumulativeSize;
    130       for (var i = 0; i < dimension; i++) centerOfMass[i] *= mult1;
    131       for (var i = 0; i < dimension; i++) centerOfMass[i] += mult2 * point[i];
     130      for(var i = 0; i < dimension; i++) centerOfMass[i] *= mult1;
     131      for(var i = 0; i < dimension; i++) centerOfMass[i] += mult2 * point[i];
    132132
    133133      // If there is space in this quad tree and it is a leaf, add the object here
    134       if (isLeaf && size < QT_NODE_CAPACITY) {
     134      if(isLeaf && size < QT_NODE_CAPACITY) {
    135135        index[size] = newIndex;
    136136        size++;
     
    140140      // Don't add duplicates for now (this is not very nice)
    141141      var anyDuplicate = false;
    142       for (uint n = 0; n < size; n++) {
     142      for(uint n = 0; n < size; n++) {
    143143        var duplicate = true;
    144         for (var d = 0; d < dimension; d++) {
    145           if (Math.Abs(point[d] - Data[index[n], d]) < double.Epsilon) continue;
     144        for(var d = 0; d < dimension; d++) {
     145          if(Math.Abs(point[d] - data[index[n], d]) < double.Epsilon) continue;
    146146          duplicate = false; break;
    147147        }
    148148        anyDuplicate = anyDuplicate | duplicate;
    149149      }
    150       if (anyDuplicate) return true;
     150      if(anyDuplicate) return true;
    151151
    152152      // Otherwise, we need to subdivide the current cell
    153       if (isLeaf) Subdivide();
     153      if(isLeaf) Subdivide();
    154154      // Find out where the point can be inserted
    155       for (var i = 0; i < noChildren; i++) {
    156         if (children[i].Insert(newIndex)) return true;
     155      for(var i = 0; i < noChildren; i++) {
     156        if(children[i].Insert(newIndex)) return true;
    157157      }
    158158
     
    165165      var newCorner = new double[dimension];
    166166      var newWidth = new double[dimension];
    167       for (var i = 0; i < noChildren; i++) {
     167      for(var i = 0; i < noChildren; i++) {
    168168        var div = 1;
    169         for (var d = 0; d < dimension; d++) {
     169        for(var d = 0; d < dimension; d++) {
    170170          newWidth[d] = .5 * boundary.GetWidth(d);
    171           if ((i / div) % 2 == 1) newCorner[d] = boundary.GetCorner(d) - .5 * boundary.GetWidth(d);
     171          if((i / div) % 2 == 1) newCorner[d] = boundary.GetCorner(d) - .5 * boundary.GetWidth(d);
    172172          else newCorner[d] = boundary.GetCorner(d) + .5 * boundary.GetWidth(d);
    173173          div *= 2;
    174174        }
    175         children[i] = new SpacePartitioningTree(this, Data, newCorner, newWidth);
     175        children[i] = new SpacePartitioningTree(this, data, newCorner, newWidth);
    176176      }
    177177
    178178      // Move existing points to correct children
    179       for (var i = 0; i < size; i++) {
     179      for(var i = 0; i < size; i++) {
    180180        var success = false;
    181         for (var j = 0; j < noChildren; j++) {
    182           if (!success) success = children[j].Insert(index[i]);
     181        for(var j = 0; j < noChildren; j++) {
     182          if(!success) success = children[j].Insert(index[i]);
    183183        }
    184184        index[i] = int.MaxValue;
     
    192192    public bool IsCorrect() {
    193193      var row = new double[dimension];
    194       for (var n = 0; n < size; n++) Buffer.BlockCopy(Data, sizeof(double) * dimension * n, row, 0, sizeof(double) * dimension);
    195       if (!boundary.ContainsPoint(row)) return false;
    196       if (isLeaf) return true;
     194      for(var n = 0; n < size; n++) Buffer.BlockCopy(data, sizeof(double) * dimension * n, row, 0, sizeof(double) * dimension);
     195      if(!boundary.ContainsPoint(row)) return false;
     196      if(isLeaf) return true;
    197197      var correct = true;
    198       for (var i = 0; i < noChildren; i++) correct = correct && children[i].IsCorrect();
     198      for(var i = 0; i < noChildren; i++) correct = correct && children[i].IsCorrect();
    199199      return correct;
    200200    }
     
    206206    public int GetAllIndices(int[] indices, int loc) {
    207207      // Gather indices in current quadrant
    208       for (var i = 0; i < size; i++) indices[loc + i] = index[i];
     208      for(var i = 0; i < size; i++) indices[loc + i] = index[i];
    209209      loc += (int)size;
    210210      // Gather indices in children
    211       if (isLeaf) return loc;
    212       for (var i = 0; i < noChildren; i++) loc = children[i].GetAllIndices(indices, loc);
     211      if(isLeaf) return loc;
     212      for(var i = 0; i < noChildren; i++) loc = children[i].GetAllIndices(indices, loc);
    213213      return loc;
    214214    }
     
    218218    }
    219219
    220     public void ComputeNonEdgeForces(int pointIndex, double theta, double[] negF, ref double sumQ)
    221     {
     220    public void ComputeNonEdgeForces(int pointIndex, double theta, double[] negF, ref double sumQ) {
    222221      // Make sure that we spend no time on empty nodes or self-interactions
    223       if (cumulativeSize == 0 || (isLeaf && size == 1 && index[0] == pointIndex)) return;
     222      if(cumulativeSize == 0 || (isLeaf && size == 1 && index[0] == pointIndex)) return;
    224223
    225224      // Compute distance between point and center-of-mass
     225      // TODO: squared distance with normalized axes is used here!
    226226      var D = .0;
    227       for (var d = 0; d < dimension; d++) buff[d] = Data[pointIndex, d] - centerOfMass[d];
    228       for (var d = 0; d < dimension; d++) D += buff[d] * buff[d];
     227      for(var d = 0; d < dimension; d++) buff[d] = data[pointIndex, d] - centerOfMass[d];
     228      for(var d = 0; d < dimension; d++) D += buff[d] * buff[d];
    229229
    230230      // Check whether we can use this node as a "summary"
    231231      var maxWidth = 0.0;
    232       for (var d = 0; d < dimension; d++) {
     232      for(var d = 0; d < dimension; d++) {
    233233        var curWidth = boundary.GetWidth(d);
    234234        maxWidth = (maxWidth > curWidth) ? maxWidth : curWidth;
    235235      }
    236       if (isLeaf || maxWidth / Math.Sqrt(D) < theta) {
     236      if(isLeaf || maxWidth / Math.Sqrt(D) < theta) {
    237237
    238238        // Compute and add t-SNE force between point and current node
     
    241241        sumQ += mult;
    242242        mult *= D;
    243         for (var d = 0; d < dimension; d++) negF[d] += mult * buff[d];
     243        for(var d = 0; d < dimension; d++) negF[d] += mult * buff[d];
    244244      } else {
    245245
    246246        // Recursively apply Barnes-Hut to children
    247         for (var i = 0; i < noChildren; i++) children[i].ComputeNonEdgeForces(pointIndex, theta, negF, ref sumQ);
     247        for(var i = 0; i < noChildren; i++) children[i].ComputeNonEdgeForces(pointIndex, theta, negF, ref sumQ);
    248248      }
    249249    }
     
    251251    public void ComputeEdgeForces(int[] rowP, int[] colP, double[] valP, int n, double[,] posF) {
    252252      // Loop over all edges in the graph
    253       for (var k = 0; k < n; k++) {
    254         for (var i = rowP[k]; i < rowP[k + 1]; i++) {
     253      for(var k = 0; k < n; k++) {
     254        for(var i = rowP[k]; i < rowP[k + 1]; i++) {
    255255
    256256          // Compute pairwise distance and Q-value
     257          // uses squared distance
    257258          var d = 1.0;
    258           for (var j = 0; j < dimension; j++) buff[j] = Data[k, j] - Data[colP[i], j];
    259           for (var j = 0; j < dimension; j++) d += buff[j] * buff[j];
     259          for(var j = 0; j < dimension; j++) buff[j] = data[k, j] - data[colP[i], j];
     260          for(var j = 0; j < dimension; j++) d += buff[j] * buff[j];
    260261          d = valP[i] / d;
    261262
    262263          // Sum positive force
    263           for (var j = 0; j < dimension; j++) posF[k, j] += d * buff[j];
     264          for(var j = 0; j < dimension; j++) posF[k, j] += d * buff[j];
    264265        }
    265266      }
     
    268269    #region Helpers
    269270    private void Fill(int n) {
    270       for (var i = 0; i < n; i++) Insert(i);
     271      for(var i = 0; i < n; i++) Insert(i);
    271272    }
    272273    private void Init(SpacePartitioningTree p, double[,] inpData, IEnumerable<double> inpCorner, IEnumerable<double> inpWidth) {
     
    274275      dimension = inpData.GetLength(1);
    275276      noChildren = 2;
    276       for (uint i = 1; i < dimension; i++) noChildren *= 2;
    277       Data = inpData;
     277      for(uint i = 1; i < dimension; i++) noChildren *= 2;
     278      data = inpData;
    278279      isLeaf = true;
    279280      size = 0;
     
    291292
    292293
    293     private class Cell  {
     294    private class Cell {
    294295      private readonly uint dimension;
    295296      private readonly double[] corner;
     
    315316      }
    316317      public bool ContainsPoint(double[] point) {
    317         for (var d = 0; d < dimension; d++)
    318           if (corner[d] - width[d] > point[d] || corner[d] + width[d] < point[d]) return false;
     318        for(var d = 0; d < dimension; d++)
     319          if(corner[d] - width[d] > point[d] || corner[d] + width[d] < point[d]) return false;
    319320        return true;
    320321      }
  • branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/TSNEStatic.cs

    r14788 r14806  
    7070    public sealed class TSNEState : DeepCloneable {
    7171      // initialized once
     72      [Storable]
    7273      public IDistance<T> distance;
     74      [Storable]
    7375      public IRandom random;
     76      [Storable]
    7477      public double perplexity;
     78      [Storable]
    7579      public bool exact;
     80      [Storable]
    7681      public int noDatapoints;
     82      [Storable]
    7783      public double finalMomentum;
     84      [Storable]
    7885      public int momSwitchIter;
     86      [Storable]
    7987      public int stopLyingIter;
     88      [Storable]
    8089      public double theta;
     90      [Storable]
    8191      public double eta;
     92      [Storable]
    8293      public int newDimensions;
    8394
    8495      // for approximate version: sparse representation of similarity/distance matrix
     96      [Storable]
    8597      public double[] valP; // similarity/distance
     98      [Storable]
    8699      public int[] rowP; // row index
     100      [Storable]
    87101      public int[] colP; // col index
    88102
    89103      // for exact version: dense representation of distance/similarity matrix
     104      [Storable]
    90105      public double[,] p;
    91106
    92107      // mapped data
     108      [Storable]
    93109      public double[,] newData;
    94110
     111      [Storable]
    95112      public int iter;
     113      [Storable]
    96114      public double currentMomentum;
    97115
    98116      // helper variables (updated in each iteration)
     117      [Storable]
    99118      public double[,] gains;
     119      [Storable]
    100120      public double[,] uY;
     121      [Storable]
    101122      public double[,] dY;
    102123
    103124      private TSNEState(TSNEState original, Cloner cloner) : base(original, cloner) {
    104       }
     125        this.distance = cloner.Clone(original.distance);
     126        this.random = cloner.Clone(original.random);
     127        this.perplexity = original.perplexity;
     128        this.exact = original.exact;
     129        this.noDatapoints = original.noDatapoints;
     130        this.finalMomentum = original.finalMomentum;
     131        this.momSwitchIter = original.momSwitchIter;
     132        this.stopLyingIter = original.stopLyingIter;
     133        this.theta = original.theta;
     134        this.eta = original.eta;
     135        this.newDimensions = original.newDimensions;
     136        if(original.valP != null) {
     137          this.valP = new double[original.valP.Length];
     138          Array.Copy(original.valP, this.valP, this.valP.Length);
     139        }
     140        if(original.rowP != null) {
     141          this.rowP = new int[original.rowP.Length];
     142          Array.Copy(original.rowP, this.rowP, this.rowP.Length);
     143        }
     144        if(original.colP != null) {
     145          this.colP = new int[original.colP.Length];
     146          Array.Copy(original.colP, this.colP, this.colP.Length);
     147        }
     148        if(original.p != null) {
     149          this.p = new double[original.p.GetLength(0), original.p.GetLength(1)];
     150          Array.Copy(original.p, this.p, this.p.Length);
     151        }
     152        this.newData = new double[original.newData.GetLength(0), original.newData.GetLength(1)];
     153        Array.Copy(original.newData, this.newData, this.newData.Length);
     154        this.iter = original.iter;
     155        this.currentMomentum = original.currentMomentum;
     156        this.gains = new double[original.gains.GetLength(0), original.gains.GetLength(1)];
     157        Array.Copy(original.gains, this.gains, this.gains.Length);
     158        this.uY = new double[original.uY.GetLength(0), original.uY.GetLength(1)];
     159        Array.Copy(original.uY, this.uY, this.uY.Length);
     160        this.dY = new double[original.dY.GetLength(0), original.dY.GetLength(1)];
     161        Array.Copy(original.dY, this.dY, this.dY.Length);
     162      }
     163
    105164      public override IDeepCloneable Clone(Cloner cloner) {
    106165        return new TSNEState(this, cloner);
     
    122181        // initialize
    123182        noDatapoints = data.Length;
    124         if (noDatapoints - 1 < 3 * perplexity) throw new ArgumentException("Perplexity too large for the number of data points!");
     183        if(noDatapoints - 1 < 3 * perplexity)
     184          throw new ArgumentException("Perplexity too large for the number of data points!");
    125185
    126186        exact = Math.Abs(theta) < double.Epsilon;
     
    129189        uY = new double[noDatapoints, newDimensions];
    130190        gains = new double[noDatapoints, newDimensions];
    131         for (var i = 0; i < noDatapoints; i++)
    132           for (var j = 0; j < newDimensions; j++)
     191        for(var i = 0; i < noDatapoints; i++)
     192          for(var j = 0; j < newDimensions; j++)
    133193            gains[i, j] = 1.0;
    134194
     
    139199
    140200        //Calculate Similarities
    141         if (exact) p = CalculateExactSimilarites(data, distance, perplexity);
     201        if(exact) p = CalculateExactSimilarites(data, distance, perplexity);
    142202        else CalculateApproximateSimilarities(data, distance, perplexity, out rowP, out colP, out valP);
    143203
    144204        // Lie about the P-values
    145         if (exact) for (var i = 0; i < noDatapoints; i++) for (var j = 0; j < noDatapoints; j++) p[i, j] *= 12.0;
    146         else for (var i = 0; i < rowP[noDatapoints]; i++) valP[i] *= 12.0;
     205        if(exact) for(var i = 0; i < noDatapoints; i++) for(var j = 0; j < noDatapoints; j++) p[i, j] *= 12.0;
     206        else for(var i = 0; i < rowP[noDatapoints]; i++) valP[i] *= 12.0;
    147207
    148208        // Initialize solution (randomly)
    149209        var rand = new NormalDistributedRandom(random, 0, 1);
    150         for (var i = 0; i < noDatapoints; i++)
    151           for (var j = 0; j < newDimensions; j++)
    152             newData[i, j] = rand.NextDouble() * .0001;  // TODO const
     210        for(var i = 0; i < noDatapoints; i++)
     211          for(var j = 0; j < newDimensions; j++)
     212            newData[i, j] = rand.NextDouble() * .0001;  // TODO const  ?
    153213      }
    154214
    155215      public double EvaluateError() {
    156         return exact ? EvaluateErrorExact(p, newData, noDatapoints, newDimensions) : EvaluateErrorApproximate(rowP, colP, valP, newData, theta);
     216        return exact ?
     217          EvaluateErrorExact(p, newData, noDatapoints, newDimensions) :
     218          EvaluateErrorApproximate(rowP, colP, valP, newData, theta);
    157219      }
    158220
     
    168230        valP = sValP;
    169231        var sumP = .0;
    170         for (var i = 0; i < rowP[data.Length]; i++) sumP += valP[i];
    171         for (var i = 0; i < rowP[data.Length]; i++) valP[i] /= sumP;
    172       }
     232        for(var i = 0; i < rowP[data.Length]; i++) sumP += valP[i];
     233        for(var i = 0; i < rowP[data.Length]; i++) valP[i] /= sumP;
     234      }
     235
    173236      private static double[,] CalculateExactSimilarites(T[] data, IDistance<T> distance, double perplexity) {
    174237        // Compute similarities
     
    176239        ComputeGaussianPerplexity(data, distance, p, perplexity);
    177240        // Symmetrize input similarities
    178         for (var n = 0; n < data.Length; n++) {
    179           for (var m = n + 1; m < data.Length; m++) {
     241        for(var n = 0; n < data.Length; n++) {
     242          for(var m = n + 1; m < data.Length; m++) {
    180243            p[n, m] += p[m, n];
    181244            p[m, n] = p[n, m];
     
    183246        }
    184247        var sumP = .0;
    185         for (var i = 0; i < data.Length; i++) for (var j = 0; j < data.Length; j++) sumP += p[i, j];
    186         for (var i = 0; i < data.Length; i++) for (var j = 0; j < data.Length; j++) p[i, j] /= sumP;
     248        for(var i = 0; i < data.Length; i++) for(var j = 0; j < data.Length; j++) sumP += p[i, j];
     249        for(var i = 0; i < data.Length; i++) for(var j = 0; j < data.Length; j++) p[i, j] /= sumP;
    187250        return p;
    188251      }
    189252
    190253      private static void ComputeGaussianPerplexity(IReadOnlyList<T> x, IDistance<T> distance, out int[] rowP, out int[] colP, out double[] valP, double perplexity, int k) {
    191         if (perplexity > k) throw new ArgumentException("Perplexity should be lower than k!");
     254        if(perplexity > k) throw new ArgumentException("Perplexity should be lower than k!");
    192255
    193256        int n = x.Count;
     
    198261        var curP = new double[n - 1];
    199262        rowP[0] = 0;
    200         for (var i = 0; i < n; i++) rowP[i + 1] = rowP[i] + k;
     263        for(var i = 0; i < n; i++) rowP[i + 1] = rowP[i] + k;
    201264
    202265        var objX = new List<IndexedItem<T>>();
    203         for (var i = 0; i < n; i++) objX.Add(new IndexedItem<T>(i, x[i]));
     266        for(var i = 0; i < n; i++) objX.Add(new IndexedItem<T>(i, x[i]));
    204267
    205268        // Build ball tree on data set
     
    207270
    208271        // Loop over all points to find nearest neighbors
    209         for (var i = 0; i < n; i++) {
     272        for(var i = 0; i < n; i++) {
    210273          IList<IndexedItem<T>> indices;
    211274          IList<double> distances;
     
    223286          // Iterate until we found a good perplexity
    224287          var iter = 0; double sumP = 0;
    225           while (!found && iter < 200) {
     288          while(!found && iter < 200) {  // TODO 200 iterations always ok?
    226289
    227290            // Compute Gaussian kernel row
    228             for (var m = 0; m < k; m++) curP[m] = Math.Exp(-beta * distances[m + 1]);
     291            for(var m = 0; m < k; m++) curP[m] = Math.Exp(-beta * distances[m + 1]); // TODO distances m+1?
    229292
    230293            // Compute entropy of current row
    231294            sumP = double.Epsilon;
    232             for (var m = 0; m < k; m++) sumP += curP[m];
     295            for(var m = 0; m < k; m++) sumP += curP[m];
    233296            var h = .0;
    234             for (var m = 0; m < k; m++) h += beta * (distances[m + 1] * curP[m]);
     297            for(var m = 0; m < k; m++) h += beta * (distances[m + 1] * curP[m]); // TODO: distances m+1?
    235298            h = h / sumP + Math.Log(sumP);
    236299
    237300            // Evaluate whether the entropy is within the tolerance level
    238301            var hdiff = h - Math.Log(perplexity);
    239             if (hdiff < tol && -hdiff < tol) {
     302            if(hdiff < tol && -hdiff < tol) {
    240303              found = true;
    241304            } else {
    242               if (hdiff > 0) {
     305              if(hdiff > 0) {
    243306                minBeta = beta;
    244                 if (maxBeta.IsAlmost(double.MaxValue) || maxBeta.IsAlmost(double.MinValue))
     307                if(maxBeta.IsAlmost(double.MaxValue) || maxBeta.IsAlmost(double.MinValue))
    245308                  beta *= 2.0;
    246309                else
     
    248311              } else {
    249312                maxBeta = beta;
    250                 if (minBeta.IsAlmost(double.MinValue) || minBeta.IsAlmost(double.MaxValue))
     313                if(minBeta.IsAlmost(double.MinValue) || minBeta.IsAlmost(double.MaxValue))
    251314                  beta /= 2.0;
    252315                else
     
    260323
    261324          // Row-normalize current row of P and store in matrix
    262           for (var m = 0; m < k; m++) curP[m] /= sumP;
    263           for (var m = 0; m < k; m++) {
     325          for(var m = 0; m < k; m++) curP[m] /= sumP;
     326          for(var m = 0; m < k; m++) {
    264327            colP[rowP[i] + m] = indices[m + 1].Index;
    265328            valP[rowP[i] + m] = curP[m];
     
    273336        int n = x.Length;
    274337        // Compute the Gaussian kernel row by row
    275         for (var i = 0; i < n; i++) {
     338        for(var i = 0; i < n; i++) {
    276339          // Initialize some variables
    277340          var found = false;
     
    284347          // Iterate until we found a good perplexity
    285348          var iter = 0;
    286           while (!found && iter < 200) {       // TODO constant
     349          while(!found && iter < 200) {       // TODO constant
    287350
    288351            // Compute Gaussian kernel row
    289             for (var m = 0; m < n; m++) p[i, m] = Math.Exp(-beta * dd[i][m]);
     352            for(var m = 0; m < n; m++) p[i, m] = Math.Exp(-beta * dd[i][m]);
    290353            p[i, i] = double.Epsilon;
    291354
    292355            // Compute entropy of current row
    293356            sumP = double.Epsilon;
    294             for (var m = 0; m < n; m++) sumP += p[i, m];
     357            for(var m = 0; m < n; m++) sumP += p[i, m];
    295358            var h = 0.0;
    296             for (var m = 0; m < n; m++) h += beta * (dd[i][m] * p[i, m]);
     359            for(var m = 0; m < n; m++) h += beta * (dd[i][m] * p[i, m]);
    297360            h = h / sumP + Math.Log(sumP);
    298361
    299362            // Evaluate whether the entropy is within the tolerance level
    300363            var hdiff = h - Math.Log(perplexity);
    301             if (hdiff < tol && -hdiff < tol) {
     364            if(hdiff < tol && -hdiff < tol) {
    302365              found = true;
    303366            } else {
    304               if (hdiff > 0) {
     367              if(hdiff > 0) {
    305368                minBeta = beta;
    306                 if (maxBeta.IsAlmost(double.MaxValue) || maxBeta.IsAlmost(double.MinValue))
     369                if(maxBeta.IsAlmost(double.MaxValue) || maxBeta.IsAlmost(double.MinValue))
    307370                  beta *= 2.0;
    308371                else
     
    310373              } else {
    311374                maxBeta = beta;
    312                 if (minBeta.IsAlmost(double.MinValue) || minBeta.IsAlmost(double.MaxValue))
     375                if(minBeta.IsAlmost(double.MinValue) || minBeta.IsAlmost(double.MaxValue))
    313376                  beta /= 2.0;
    314377                else
     
    322385
    323386          // Row normalize P
    324           for (var m = 0; m < n; m++) p[i, m] /= sumP;
     387          for(var m = 0; m < n; m++) p[i, m] /= sumP;
    325388        }
    326389      }
    327390
    328391      private static double[][] ComputeDistances(T[] x, IDistance<T> distance) {
    329         return x.Select(m => x.Select(n => distance.Get(m, n)).ToArray()).ToArray();
    330       }
    331 
    332 
     392        var res = new double[x.Length][];
     393        for(int r = 0; r < x.Length; r++) {
     394          var rowV = new double[x.Length];
     395          // all distances must be symmetric
     396          for(int c = 0; c < r; c++) {
     397            rowV[c] = res[c][r];
     398          }
     399          rowV[r] = 0.0; // distance to self is zero for all distances
     400          for(int c = r + 1; c < x.Length; c++) {
     401            rowV[c] = distance.Get(x[r], x[c]);
     402          }
     403          res[r] = rowV;
     404        }
     405        return res;
     406        // return x.Select(m => x.Select(n => distance.Get(m, n)).ToArray()).ToArray();
     407      }
    333408
    334409      private static double EvaluateErrorExact(double[,] p, double[,] y, int n, int d) {
     
    336411        var dd = new double[n, n];
    337412        var q = new double[n, n];
    338         ComputeSquaredEuclideanDistance(y, n, d, dd);
     413        ComputeSquaredEuclideanDistance(y, n, d, dd); // TODO: we use Euclidian distance regardless of the actual distance function
    339414
    340415        // Compute Q-matrix and normalization sum
    341416        var sumQ = double.Epsilon;
    342         for (var n1 = 0; n1 < n; n1++) {
    343           for (var m = 0; m < n; m++) {
    344             if (n1 != m) {
     417        for(var n1 = 0; n1 < n; n1++) {
     418          for(var m = 0; m < n; m++) {
     419            if(n1 != m) {
    345420              q[n1, m] = 1 / (1 + dd[n1, m]);
    346421              sumQ += q[n1, m];
     
    348423          }
    349424        }
    350         for (var i = 0; i < n; i++) for (var j = 0; j < n; j++) q[i, j] /= sumQ;
     425        for(var i = 0; i < n; i++) for(var j = 0; j < n; j++) q[i, j] /= sumQ;
    351426
    352427        // Sum t-SNE error
    353428        var c = .0;
    354         for (var i = 0; i < n; i++)
    355           for (var j = 0; j < n; j++) {
     429        for(var i = 0; i < n; i++)
     430          for(var j = 0; j < n; j++) {
    356431            c += p[i, j] * Math.Log((p[i, j] + float.Epsilon) / (q[i, j] + float.Epsilon));
    357432          }
    358433        return c;
    359434      }
     435
     436      // TODO: there seems to be a bug in the error approximation.
     437      // The mapping of the approximate tSNE looks good but the error curve never changes.
    360438      private static double EvaluateErrorApproximate(IReadOnlyList<int> rowP, IReadOnlyList<int> colP, IReadOnlyList<double> valP, double[,] y, double theta) {
    361439        // Get estimate of normalization term
     
    365443        var buff = new double[d];
    366444        double sumQ = 0.0;
    367         for (var i = 0; i < n; i++) tree.ComputeNonEdgeForces(i, theta, buff, ref sumQ);
     445        for(var i = 0; i < n; i++) tree.ComputeNonEdgeForces(i, theta, buff, ref sumQ);
    368446
    369447        // Loop over all edges to compute t-SNE error
    370448        var c = .0;
    371         for (var k = 0; k < n; k++) {
    372           for (var i = rowP[k]; i < rowP[k + 1]; i++) {
     449        for(var k = 0; k < n; k++) {
     450          for(var i = rowP[k]; i < rowP[k + 1]; i++) {
    373451            var q = .0;
    374             for (var j = 0; j < d; j++) buff[j] = y[k, j];
    375             for (var j = 0; j < d; j++) buff[j] -= y[colP[i], j];
    376             for (var j = 0; j < d; j++) q += buff[j] * buff[j];
     452            for(var j = 0; j < d; j++) buff[j] = y[k, j];
     453            for(var j = 0; j < d; j++) buff[j] -= y[colP[i], j];
     454            for(var j = 0; j < d; j++) q += buff[j] * buff[j];     // TODO: squared error is used here!
    377455            q = 1.0 / (1.0 + q) / sumQ;
    378456            c += valP[i] * Math.Log((valP[i] + float.Epsilon) / (q + float.Epsilon));
     
    386464        var n = rowP.Count - 1;
    387465        var rowCounts = new int[n];
    388         for (var j = 0; j < n; j++) {
    389           for (var i = rowP[j]; i < rowP[j + 1]; i++) {
     466        for(var j = 0; j < n; j++) {
     467          for(var i = rowP[j]; i < rowP[j + 1]; i++) {
    390468
    391469            // Check whether element (col_P[i], n) is present
    392470            var present = false;
    393             for (var m = rowP[colP[i]]; m < rowP[colP[i] + 1]; m++) {
    394               if (colP[m] == j) present = true;
     471            for(var m = rowP[colP[i]]; m < rowP[colP[i] + 1]; m++) {
     472              if(colP[m] == j) present = true;
    395473            }
    396             if (present) rowCounts[j]++;
     474            if(present) rowCounts[j]++;
    397475            else {
    398476              rowCounts[j]++;
     
    402480        }
    403481        var noElem = 0;
    404         for (var i = 0; i < n; i++) noElem += rowCounts[i];
     482        for(var i = 0; i < n; i++) noElem += rowCounts[i];
    405483
    406484        // Allocate memory for symmetrized matrix
     
    411489        // Construct new row indices for symmetric matrix
    412490        symRowP[0] = 0;
    413         for (var i = 0; i < n; i++) symRowP[i + 1] = symRowP[i] + rowCounts[i];
     491        for(var i = 0; i < n; i++) symRowP[i + 1] = symRowP[i] + rowCounts[i];
    414492
    415493        // Fill the result matrix
    416494        var offset = new int[n];
    417         for (var j = 0; j < n; j++) {
    418           for (var i = rowP[j]; i < rowP[j + 1]; i++) {                                  // considering element(n, colP[i])
     495        for(var j = 0; j < n; j++) {
     496          for(var i = rowP[j]; i < rowP[j + 1]; i++) {                                  // considering element(n, colP[i])
    419497
    420498            // Check whether element (col_P[i], n) is present
    421499            var present = false;
    422             for (var m = rowP[colP[i]]; m < rowP[colP[i] + 1]; m++) {
    423               if (colP[m] != j) continue;
     500            for(var m = rowP[colP[i]]; m < rowP[colP[i] + 1]; m++) {
     501              if(colP[m] != j) continue;
    424502              present = true;
    425               if (j > colP[i]) continue; // make sure we do not add elements twice
     503              if(j > colP[i]) continue; // make sure we do not add elements twice
    426504              symColP[symRowP[j] + offset[j]] = colP[i];
    427505              symColP[symRowP[colP[i]] + offset[colP[i]]] = j;
     
    431509
    432510            // If (colP[i], n) is not present, there is no addition involved
    433             if (!present) {
     511            if(!present) {
    434512              symColP[symRowP[j] + offset[j]] = colP[i];
    435513              symColP[symRowP[colP[i]] + offset[colP[i]]] = j;
     
    439517
    440518            // Update offsets
    441             if (present && (j > colP[i])) continue;
     519            if(present && (j > colP[i])) continue;
    442520            offset[j]++;
    443             if (colP[i] != j) offset[colP[i]]++;
    444           }
    445         }
    446 
    447         // Divide the result by two
    448         for (var i = 0; i < noElem; i++) symValP[i] /= 2.0;
     521            if(colP[i] != j) offset[colP[i]]++;
     522          }
     523        }
     524
     525        for(var i = 0; i < noElem; i++) symValP[i] /= 2.0;
    449526      }
    450527
     
    459536
    460537    public static double[,] Iterate(TSNEState state) {
    461       if (state.exact)
     538      if(state.exact)
    462539        ComputeExactGradient(state.p, state.newData, state.noDatapoints, state.newDimensions, state.dY);
    463540      else
     
    465542
    466543      // Update gains
    467       for (var i = 0; i < state.noDatapoints; i++) {
    468         for (var j = 0; j < state.newDimensions; j++) {
     544      for(var i = 0; i < state.noDatapoints; i++) {
     545        for(var j = 0; j < state.newDimensions; j++) {
    469546          state.gains[i, j] = Math.Sign(state.dY[i, j]) != Math.Sign(state.uY[i, j])
    470547            ? state.gains[i, j] + .2
    471548            : state.gains[i, j] * .8; // 20% up or 20% down // TODO: +0.2?!
    472549
    473           if (state.gains[i, j] < .01) state.gains[i, j] = .01;
     550          if(state.gains[i, j] < .01) state.gains[i, j] = .01; // TODO why limit the gains?
    474551        }
    475552      }
     
    477554
    478555      // Perform gradient update (with momentum and gains)
    479       for (var i = 0; i < state.noDatapoints; i++)
    480         for (var j = 0; j < state.newDimensions; j++)
     556      for(var i = 0; i < state.noDatapoints; i++)
     557        for(var j = 0; j < state.newDimensions; j++)
    481558          state.uY[i, j] = state.currentMomentum * state.uY[i, j] - state.eta * state.gains[i, j] * state.dY[i, j];
    482559
    483       for (var i = 0; i < state.noDatapoints; i++)
    484         for (var j = 0; j < state.newDimensions; j++)
     560      for(var i = 0; i < state.noDatapoints; i++)
     561        for(var j = 0; j < state.newDimensions; j++)
    485562          state.newData[i, j] = state.newData[i, j] + state.uY[i, j];
    486563
     
    489566      // Stop lying about the P-values after a while, and switch momentum
    490567
    491       if (state.iter == state.stopLyingIter) {
    492         if (state.exact)
    493           for (var i = 0; i < state.noDatapoints; i++) for (var j = 0; j < state.noDatapoints; j++) state.p[i, j] /= 12.0;                                   //XXX why 12?
     568      if(state.iter == state.stopLyingIter) {
     569        if(state.exact)
     570          for(var i = 0; i < state.noDatapoints; i++) for(var j = 0; j < state.noDatapoints; j++) state.p[i, j] /= 12.0;                                   //XXX why 12?
    494571        else
    495           for (var i = 0; i < state.rowP[state.noDatapoints]; i++) state.valP[i] /= 12.0;                       // XXX are we not scaling all values?
    496       }
    497 
    498       if (state.iter == state.momSwitchIter)
     572          for(var i = 0; i < state.rowP[state.noDatapoints]; i++) state.valP[i] /= 12.0;                       // XXX are we not scaling all values?
     573      }
     574
     575      if(state.iter == state.momSwitchIter)
    499576        state.currentMomentum = state.finalMomentum;
    500577
     
    511588      tree.ComputeEdgeForces(rowP, colP, valP, n, posF);
    512589      var row = new double[d];
    513       for (var n1 = 0; n1 < n; n1++) {
     590      for(var n1 = 0; n1 < n; n1++) {
    514591        Buffer.BlockCopy(negF, (sizeof(double) * n1 * d), row, 0, d);
    515592        tree.ComputeNonEdgeForces(n1, theta, row, ref sumQ);
     
    517594
    518595      // Compute final t-SNE gradient
    519       for (var i = 0; i < n; i++)
    520         for (var j = 0; j < d; j++) {
     596      for(var i = 0; i < n; i++)
     597        for(var j = 0; j < d; j++) {
    521598          dC[i, j] = posF[i, j] - negF[i, j] / sumQ;
    522599        }
     
    526603
    527604      // Make sure the current gradient contains zeros
    528       for (var i = 0; i < n; i++) for (var j = 0; j < d; j++) dC[i, j] = 0.0;
     605      for(var i = 0; i < n; i++) for(var j = 0; j < d; j++) dC[i, j] = 0.0;
    529606
    530607      // Compute the squared Euclidean distance matrix
    531608      var dd = new double[n, n];
    532       ComputeSquaredEuclideanDistance(y, n, d, dd);
     609      ComputeSquaredEuclideanDistance(y, n, d, dd); // TODO: we use Euclidian distance regardless which distance function is actually set!
    533610
    534611      // Compute Q-matrix and normalization sum
    535612      var q = new double[n, n];
    536613      var sumQ = .0;
    537       for (var n1 = 0; n1 < n; n1++) {
    538         for (var m = 0; m < n; m++) {
    539           if (n1 == m) continue;
     614      for(var n1 = 0; n1 < n; n1++) {
     615        for(var m = 0; m < n; m++) {
     616          if(n1 == m) continue;
    540617          q[n1, m] = 1 / (1 + dd[n1, m]);
    541618          sumQ += q[n1, m];
     
    544621
    545622      // Perform the computation of the gradient
    546       for (var n1 = 0; n1 < n; n1++) {
    547         for (var m = 0; m < n; m++) {
    548           if (n1 == m) continue;
     623      for(var n1 = 0; n1 < n; n1++) {
     624        for(var m = 0; m < n; m++) {
     625          if(n1 == m) continue;
    549626          var mult = (p[n1, m] - q[n1, m] / sumQ) * q[n1, m];
    550           for (var d1 = 0; d1 < d; d1++) {
     627          for(var d1 = 0; d1 < d; d1++) {
    551628            dC[n1, d1] += (y[n1, d1] - y[m, d1]) * mult;
    552629          }
     
    557634    private static void ComputeSquaredEuclideanDistance(double[,] x, int n, int d, double[,] dd) {
    558635      var dataSums = new double[n];
    559       for (var i = 0; i < n; i++) {
    560         for (var j = 0; j < d; j++) {
     636      for(var i = 0; i < n; i++) {
     637        for(var j = 0; j < d; j++) {
    561638          dataSums[i] += x[i, j] * x[i, j];
    562639        }
    563640      }
    564       for (var i = 0; i < n; i++) {
    565         for (var m = 0; m < n; m++) {
     641      for(var i = 0; i < n; i++) {
     642        for(var m = 0; m < n; m++) {
    566643          dd[i, m] = dataSums[i] + dataSums[m];
    567644        }
    568645      }
    569       for (var i = 0; i < n; i++) {
     646      for(var i = 0; i < n; i++) {
    570647        dd[i, i] = 0.0;
    571         for (var m = i + 1; m < n; m++) {
     648        for(var m = i + 1; m < n; m++) {
    572649          dd[i, m] = 0.0;
    573           for (var j = 0; j < d; j++) {
     650          for(var j = 0; j < d; j++) {
    574651            dd[i, m] += (x[i, j] - x[m, j]) * (x[i, j] - x[m, j]);
    575652          }
     
    584661      var d = x.GetLength(1);
    585662      var mean = new double[d];
    586       for (var i = 0; i < n; i++) {
    587         for (var j = 0; j < d; j++) {
     663      for(var i = 0; i < n; i++) {
     664        for(var j = 0; j < d; j++) {
    588665          mean[j] += x[i, j];
    589666        }
    590667      }
    591       for (var i = 0; i < d; i++) {
     668      for(var i = 0; i < d; i++) {
    592669        mean[i] /= n;
    593670      }
    594671      // Subtract data mean
    595       for (var i = 0; i < n; i++) {
    596         for (var j = 0; j < d; j++) {
     672      for(var i = 0; i < n; i++) {
     673        for(var j = 0; j < d; j++) {
    597674          x[i, j] -= mean[j];
    598675        }
Note: See TracChangeset for help on using the changeset viewer.