Changeset 15207 for trunk/sources/HeuristicLab.Algorithms.DataAnalysis

Timestamp:

07/12/17 15:17:41 (7 years ago)

Author:

bwerth

Message:

#2700 worked on TSNE (mostly removing unused code)

Location:

trunk/sources/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE

Files:

• Distances/CosineDistance.cs (moved) (moved from trunk/sources/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/Distances/InnerProductDistance.cs) (1 diff)
• Distances/DistanceBase.cs (modified) (2 diffs)
• Distances/EuclideanDistance.cs (modified) (2 diffs)
• Distances/IDistance.cs (added)
• Distances/ManhattanDistance.cs (modified) (2 diffs)
• Interfaces (deleted)
• PriorityQueue.cs (modified) (8 diffs)
• SpacePartitioningTree.cs (modified) (14 diffs)
• TSNEAlgorithm.cs (modified) (3 diffs)
• TSNEStatic.cs (modified) (40 diffs)
• VantagePointTree.cs (modified) (3 diffs)

trunk/sources/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/Distances/CosineDistance.cs

@@ -34 +34 @@
   /// </summary>
   [StorableClass]
-  [Item("InnerProductDistance", "The angluar distance as defined as a normalized distance measure dependent on the angle between two vectors.\nIt is designed for vectors with all positive coordinates")]
-  public class InnerProductDistance : DistanceBase<IEnumerable<double>> {
+  [Item("CosineDistance", "The angluar distance as defined as a normalized distance measure dependent on the angle between two vectors.\nIt is designed for vectors with all positive coordinates")]
+  public class CosineDistance : DistanceBase<IEnumerable<double>> {
 
-    #region HLConstructors
+    #region HLConstructors & Cloning
     [StorableConstructor]
-    protected InnerProductDistance(bool deserializing) : base(deserializing) { }
-    protected InnerProductDistance(InnerProductDistance original, Cloner cloner)
+    protected CosineDistance(bool deserializing) : base(deserializing) { }
+    protected CosineDistance(CosineDistance original, Cloner cloner)
       : base(original, cloner) { }
-    public InnerProductDistance() { }
+    public CosineDistance() { }
     public override IDeepCloneable Clone(Cloner cloner) {
-      return new InnerProductDistance(this, cloner);
+      return new CosineDistance(this, cloner);
     }
     #endregion
 
     #region statics
-    public static double GetDistance(IEnumerable<double> point1, IEnumerable<double> point2) {
-      var xs = point1.GetEnumerator();
-      var ys = point2.GetEnumerator();
-      var sum = 0.0;
-      while(xs.MoveNext() & ys.MoveNext()) {
-        if(xs.Current < 0 || ys.Current < 0) throw new ArgumentException("Inner product distance is only defined for vectors with non-negative elements");
-        sum += xs.Current * ys.Current;
+    public static double GetDistance(IReadOnlyList<double> point1, IReadOnlyList<double> point2) {
+      if (point1.Count != point2.Count) throw new ArgumentException("Cosine distance not defined on vectors of different length");
+      var innerprod = 0.0;
+      var length1 = 0.0;
+      var length2 = 0.0;
+
+      for (var i = 0; i < point1.Count; i++) {
+        double d1 = point1[i], d2 = point2[i];
+        innerprod += d1 * d2;
+        length1 += d1 * d1;
+        length2 += d2 * d2;
       }
-      if(xs.MoveNext() || ys.MoveNext()) throw new ArgumentException("Enumerables contain a different number of elements");
-      return sum;
+      var l = Math.Sqrt(length1 * length2);
+      if (l.IsAlmost(0)) throw new ArgumentException("Cosine distance is not defined on vectors of length 0");
+      return 1 - innerprod / l;
     }
     #endregion
     public override double Get(IEnumerable<double> a, IEnumerable<double> b) {
-      return GetDistance(a, b);
+      return GetDistance(a.ToArray(), b.ToArray());
     }
   }

trunk/sources/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/Distances/DistanceBase.cs

@@ -30 +30 @@
   public abstract class DistanceBase<T> : Item, IDistance<T> {
 
-    #region HLConstructors & Boilerplate
+    #region HLConstructors & Cloning
     [StorableConstructor]
     protected DistanceBase(bool deserializing) : base(deserializing) { }
@@ -42 +42 @@
       return new DistanceComparer(item, this);
     }
-
 
     public double Get(object x, object y) {

trunk/sources/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/Distances/EuclideanDistance.cs

@@ -32 +32 @@
   public class EuclideanDistance : DistanceBase<IEnumerable<double>> {
 
-    #region HLConstructors & Boilerplate
+    #region HLConstructors & Cloning
     [StorableConstructor]
     protected EuclideanDistance(bool deserializing) : base(deserializing) { }
@@ -40 +40 @@
     #endregion
 
-    public static double GetDistance(double[] point1, double[] point2) {
-      if (point1.Length != point2.Length) throw new ArgumentException("Euclidean distance not defined on vectors of different length");
-      return Math.Sqrt(point1.Zip(point2, (a1, b1) => (a1 - b1) * (a1 - b1)).Sum());
+    public static double GetDistance(IReadOnlyList<double> point1, IReadOnlyList<double> point2) {
+      if (point1.Count != point2.Count) throw new ArgumentException("Euclidean distance not defined on vectors of different length");
+      var sum = 0.0;
+      for (var i = 0; i < point1.Count; i++) {
+        var d = point1[i] - point2[i];
+        sum += d * d;
+      }
+
+      return Math.Sqrt(sum);
     }
 

trunk/sources/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/Distances/ManhattanDistance.cs

@@ -32 +32 @@
   public class ManhattanDistance : DistanceBase<IEnumerable<double>> {
 
-    #region HLConstructors & Boilerplate
+    #region HLConstructors & Cloning
     [StorableConstructor]
     protected ManhattanDistance(bool deserializing) : base(deserializing) { }
@@ -47 +47 @@
     public static double GetDistance(double[] point1, double[] point2) {
       if (point1.Length != point2.Length) throw new ArgumentException("Manhattan distance not defined on vectors of different length");
-      return point1.Zip(point2, (a1, b1) => Math.Abs(a1 - b1)).Sum();
+      var sum = 0.0;
+      for (var i = 0; i < point1.Length; i++)
+        sum += Math.Abs(point1[i] + point2[i]);
+      return sum;
     }
 

trunk/sources/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/PriorityQueue.cs

@@ -25 +25 @@
 
 namespace HeuristicLab.Algorithms.DataAnalysis {
-  // TODO: move to HeuristicLab.Collections
 
   // Code Taken from branch: RoutePlanning, Heap4
@@ -63 +62 @@
 
     public KeyValuePair<TK, TV> PeekMin() {
-      if(size == 0) {
-        throw new InvalidOperationException("Heap is empty");
-      }
+      if (size == 0) throw new InvalidOperationException("Heap is empty");
       return new KeyValuePair<TK, TV>(data[1].Key, data[1].Value);
     }
@@ -86 +83 @@
       size = sz - 1;
       finalLayerDist++;
-      if(finalLayerDist == finalLayerSize) {
+      if (finalLayerDist == finalLayerSize) {
         finalLayerSize >>= 2;
         finalLayerDist = 0;
       }
 
-      while(succ < sz) {
+      while (succ < sz) {
         var minKey = data[succ].Key;
         var delta = 0;
 
-        var otherKey = data[succ + 1].Key;
-        if(otherKey.CompareTo(minKey) < 0) {
+        for (var i = 1; i <= 3; i++) {
+          var otherKey = data[succ + i].Key;
+          if (otherKey.CompareTo(minKey) >= 0) continue;
           minKey = otherKey;
-          delta = 1;
+          delta = i;
         }
-        otherKey = data[succ + 2].Key;
-        if(otherKey.CompareTo(minKey) < 0) {
-          minKey = otherKey;
-          delta = 2;
-        }
-        otherKey = data[succ + 3].Key;
-        if(otherKey.CompareTo(minKey) < 0) {
-          minKey = otherKey;
-          delta = 3;
-        }
+        
         succ += delta;
         layerPos += delta;
@@ -136 +125 @@
       pred = pred - layerDist; // finally preds index
 
-      while(data[pred].Key.CompareTo(bubble) > 0) {  // must terminate since inf at root
+      while (data[pred].Key.CompareTo(bubble) > 0) {  // must terminate since inf at root
         data[hole].Key = data[pred].Key;
         data[hole].Value = data[pred].Value;
@@ -149 +138 @@
       data[hole].Key = bubble;
       data[hole].Value = data[sz].Value;
-
       data[sz].Key = Supremum; // mark as deleted
     }
@@ -158 +146 @@
       finalLayerDist--;
 
-      if(finalLayerDist == -1) { // layer full
-                                 // start next layer
+      if (finalLayerDist == -1) { // layer full
+                                  // start next layer
         finalLayerSize <<= 2;
         finalLayerDist = finalLayerSize - 1;
@@ -172 +160 @@
       pred = pred - layerDist; // finally preds index
       var predKey = data[pred].Key;
-      while(predKey.CompareTo(key) > 0) {
+      while (predKey.CompareTo(key) > 0) {
         data[hole].Key = predKey;
         data[hole].Value = data[pred].Value;
@@ -195 +183 @@
       var sup = Supremum;
       var cap = capacity;
-      for(var i = 1; i <= cap; i++) {
+      for (var i = 1; i <= cap; i++) {
         data[i].Key = sup;
       }

trunk/sources/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/SpacePartitioningTree.cs

@@ -56 +56 @@
 using System;
 using System.Collections.Generic;
-using System.Linq;
-using HeuristicLab.Common;
 
 namespace HeuristicLab.Algorithms.DataAnalysis {
@@ -63 +61 @@
   /// Space partitioning tree (SPTree)
   /// </summary>
-  public class SpacePartitioningTree : ISpacePartitioningTree {
-    private const uint QT_NODE_CAPACITY = 1;
-
-    private double[] buff;
-    private SpacePartitioningTree parent;
+  internal class SpacePartitioningTree {
+    private const uint QtNodeCapacity = 1;
+
+    #region Fields
     private int dimension;
     private bool isLeaf;
@@ -80 +77 @@
     // Indices in this space-partitioning tree node, corresponding center-of-mass, and list of all children
     private double[] centerOfMass;
-    private readonly int[] index = new int[QT_NODE_CAPACITY];
+    private readonly int[] index = new int[QtNodeCapacity];
 
     // Children
     private SpacePartitioningTree[] children;
     private uint noChildren;
+    #endregion
 
     public SpacePartitioningTree(double[,] inpData) {
@@ -104 +102 @@
       var width = new double[d];
       for (var i = 0; i < d; i++) width[i] = Math.Max(maxY[i] - meanY[i], meanY[i] - minY[i]) + 1e-5;
-      Init(null, inpData, meanY, width);
+      Init(inpData, meanY, width);
       Fill(n);
     }
 
-    public SpacePartitioningTree(double[,] inpData, IEnumerable<double> impCorner, IEnumerable<double> impWith) {
-      Init(null, inpData, impCorner, impWith);
-    }
-    public SpacePartitioningTree(SpacePartitioningTree parent, double[,] inpData, IEnumerable<double> impCorner, IEnumerable<double> impWith) {
-      Init(parent, inpData, impCorner, impWith);
-    }
-
-    public ISpacePartitioningTree GetParent() {
-      return parent;
+    private SpacePartitioningTree(double[,] inpData, IEnumerable<double> impCorner, IEnumerable<double> impWith) {
+      Init(inpData, impCorner, impWith);
     }
 
@@ -132 +123 @@
 
       // If there is space in this quad tree and it is a leaf, add the object here
-      if (isLeaf && size < QT_NODE_CAPACITY) {
+      if (isLeaf && size < QtNodeCapacity) {
         index[size] = newIndex;
         size++;
@@ -138 +129 @@
       }
 
-      // Don't add duplicates for now (this is not very nice)
+      // Don't add duplicates
       var anyDuplicate = false;
       for (uint n = 0; n < size; n++) {
@@ -161 +152 @@
     }
 
-    public void Subdivide() {
-      // Create new children
-      var newCorner = new double[dimension];
-      var newWidth = new double[dimension];
-      for (var i = 0; i < noChildren; i++) {
-        var div = 1;
-        for (var d = 0; d < dimension; d++) {
-          newWidth[d] = .5 * boundary.GetWidth(d);
-          if ((i / div) % 2 == 1) newCorner[d] = boundary.GetCorner(d) - .5 * boundary.GetWidth(d);
-          else newCorner[d] = boundary.GetCorner(d) + .5 * boundary.GetWidth(d);
-          div *= 2;
-        }
-        children[i] = new SpacePartitioningTree(this, data, newCorner, newWidth);
-      }
-
-      // Move existing points to correct children
-      for (var i = 0; i < size; i++) {
-        var success = false;
-        for (var j = 0; j < noChildren; j++) {
-          if (!success) success = children[j].Insert(index[i]);
-        }
-        index[i] = -1; // as in tSNE implementation by van der Maaten
-      }
-
-      // Empty parent node
-      size = 0;
-      isLeaf = false;
-    }
-
-    public bool IsCorrect() {
-      var row = new double[dimension];
-      for (var n = 0; n < size; n++) {
-        Buffer.BlockCopy(data, sizeof(double) * dimension * index[n], row, 0, sizeof(double) * dimension);
-        if (!boundary.ContainsPoint(row)) return false;
-      }
-      if (isLeaf) return true;
-      var correct = true;
-      for (var i = 0; i < noChildren; i++) correct = correct && children[i].IsCorrect();
-      return correct;
-    }
-
-    public void GetAllIndices(int[] indices) {
-      GetAllIndices(indices, 0);
-    }
-
-    public int GetAllIndices(int[] indices, int loc) {
-      // Gather indices in current quadrant
-      for (var i = 0; i < size; i++) indices[loc + i] = index[i];
-      loc += (int)size;
-      // Gather indices in children
-      if (isLeaf) return loc;
-      for (var i = 0; i < noChildren; i++) loc = children[i].GetAllIndices(indices, loc);
-      return loc;
-    }
-
-    public int GetDepth() {
-      return isLeaf ? 1 : 1 + children.Max(x => x.GetDepth());
-    }
-
     public void ComputeNonEdgeForces(int pointIndex, double theta, double[] negF, ref double sumQ) {
       // Make sure that we spend no time on empty nodes or self-interactions
@@ -226 +158 @@
       // Compute distance between point and center-of-mass
       var D = .0;
+      var buff = new double[dimension];
       for (var d = 0; d < dimension; d++) buff[d] = data[pointIndex, d] - centerOfMass[d];
       for (var d = 0; d < dimension; d++) D += buff[d] * buff[d];
@@ -233 +166 @@
       for (var d = 0; d < dimension; d++) {
         var curWidth = boundary.GetWidth(d);
-        maxWidth = (maxWidth > curWidth) ? maxWidth : curWidth;
+        maxWidth = maxWidth > curWidth ? maxWidth : curWidth;
       }
       if (isLeaf || maxWidth / Math.Sqrt(D) < theta) {
@@ -250 +183 @@
     }
 
-    // does not use the tree
-    public void ComputeEdgeForces(int[] rowP, int[] colP, double[] valP, int n, double[,] posF) {
+    public static void ComputeEdgeForces(int[] rowP, int[] colP, double[] valP, int n, double[,] posF, double[,] data, int dimension) {
       // Loop over all edges in the graph
       for (var k = 0; k < n; k++) {
@@ -259 +191 @@
           // uses squared distance
           var d = 1.0;
+          var buff = new double[dimension];
           for (var j = 0; j < dimension; j++) buff[j] = data[k, j] - data[colP[i], j];
           for (var j = 0; j < dimension; j++) d += buff[j] * buff[j];
@@ -273 +206 @@
       for (var i = 0; i < n; i++) Insert(i);
     }
-    private void Init(SpacePartitioningTree p, double[,] inpData, IEnumerable<double> inpCorner, IEnumerable<double> inpWidth) {
-      parent = p;
+
+    private void Init(double[,] inpData, IEnumerable<double> inpCorner, IEnumerable<double> inpWidth) {
       dimension = inpData.GetLength(1);
       noChildren = 2;
@@ -283 +216 @@
       cumulativeSize = 0;
       boundary = new Cell((uint)dimension);
+
       inpCorner.ForEach((i, x) => boundary.SetCorner(i, x));
       inpWidth.ForEach((i, x) => boundary.SetWidth(i, x));
@@ -288 +222 @@
       children = new SpacePartitioningTree[noChildren];
       centerOfMass = new double[dimension];
-      buff = new double[dimension];
-
+    }
+
+    private void Subdivide() {
+      // Create new children
+      var newCorner = new double[dimension];
+      var newWidth = new double[dimension];
+      for (var i = 0; i < noChildren; i++) {
+        var div = 1;
+        for (var d = 0; d < dimension; d++) {
+          newWidth[d] = .5 * boundary.GetWidth(d);
+          if (i / div % 2 == 1) newCorner[d] = boundary.GetCorner(d) - .5 * boundary.GetWidth(d);
+          else newCorner[d] = boundary.GetCorner(d) + .5 * boundary.GetWidth(d);
+          div *= 2;
+        }
+        children[i] = new SpacePartitioningTree(data, newCorner, newWidth);
+      }
+
+      // Move existing points to correct children
+      for (var i = 0; i < size; i++) {
+        var success = false;
+        for (var j = 0; j < noChildren; j++) {
+          if (!success) success = children[j].Insert(index[i]);
+        }
+        index[i] = -1; // as in tSNE implementation by van der Maaten
+      }
+      // Empty parent node
+      size = 0;
+      isLeaf = false;
     }
     #endregion
-
 
     private class Cell {

trunk/sources/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/TSNEAlgorithm.cs

@@ -41 +41 @@
   /// </summary>
   [Item("tSNE", "t-distributed stochastic neighbourhood embedding projects the data in a low " +
-                "dimensional space to allow visual cluster identification.")]
+                "dimensional space to allow visual cluster identification. Implemented similar to: https://lvdmaaten.github.io/tsne/#implementations (Barnes-Hut t-SNE). Described in : https://lvdmaaten.github.io/publications/papers/JMLR_2014.pdf")]
   [Creatable(CreatableAttribute.Categories.DataAnalysis, Priority = 100)]
   [StorableClass]
@@ -198 +198 @@
 
     private TSNEAlgorithm(TSNEAlgorithm original, Cloner cloner) : base(original, cloner) {
-      if(original.dataRowNames!=null)
-      this.dataRowNames = new Dictionary<string, List<int>>(original.dataRowNames);
+      if (original.dataRowNames != null)
+        this.dataRowNames = new Dictionary<string, List<int>>(original.dataRowNames);
       if (original.dataRows != null)
         this.dataRows = original.dataRows.ToDictionary(kvp => kvp.Key, kvp => cloner.Clone(kvp.Value));
@@ -226 +226 @@
       Parameters.Add(new FixedValueParameter<BoolValue>(SetSeedRandomlyParameterName, "If the seed should be random.", new BoolValue(true)));
       Parameters.Add(new FixedValueParameter<IntValue>(SeedParameterName, "The seed used if it should not be random.", new IntValue(0)));
-      Parameters.Add(new FixedValueParameter<StringValue>(ClassesParameterName, "name of the column specifying the class lables of each data point. If the label column can not be found training/test is used as labels.", new StringValue("none")));
+      Parameters.Add(new FixedValueParameter<StringValue>(ClassesParameterName, "Name of the column specifying the class lables of each data point. If the label column can not be found training/test is used as labels.", new StringValue("none")));
       Parameters.Add(new FixedValueParameter<BoolValue>(NormalizationParameterName, "Whether the data should be zero centered and have variance of 1 for each variable, so different scalings are ignored.", new BoolValue(true)));
       Parameters.Add(new FixedValueParameter<IntValue>(UpdateIntervalParameterName, "", new IntValue(50)));

trunk/sources/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/TSNEStatic.cs

@@ -56 +56 @@
 using System;
 using System.Collections.Generic;
-using System.Linq;
 using HeuristicLab.Collections;
 using HeuristicLab.Common;
 using HeuristicLab.Core;
-using HeuristicLab.Optimization;
 using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
 using HeuristicLab.Random;
@@ -70 +68 @@
     [StorableClass]
     public sealed class TSNEState : DeepCloneable {
+      #region Storables
       // initialized once
       [Storable]
@@ -122 +121 @@
       [Storable]
       public double[,] dY;
-
+      #endregion
+
+      #region Constructors & Cloning
       private TSNEState(TSNEState original, Cloner cloner) : base(original, cloner) {
-        this.distance = cloner.Clone(original.distance);
-        this.random = cloner.Clone(original.random);
-        this.perplexity = original.perplexity;
-        this.exact = original.exact;
-        this.noDatapoints = original.noDatapoints;
-        this.finalMomentum = original.finalMomentum;
-        this.momSwitchIter = original.momSwitchIter;
-        this.stopLyingIter = original.stopLyingIter;
-        this.theta = original.theta;
-        this.eta = original.eta;
-        this.newDimensions = original.newDimensions;
-        if(original.valP != null) {
-          this.valP = new double[original.valP.Length];
-          Array.Copy(original.valP, this.valP, this.valP.Length);
-        }
-        if(original.rowP != null) {
-          this.rowP = new int[original.rowP.Length];
-          Array.Copy(original.rowP, this.rowP, this.rowP.Length);
-        }
-        if(original.colP != null) {
-          this.colP = new int[original.colP.Length];
-          Array.Copy(original.colP, this.colP, this.colP.Length);
-        }
-        if(original.p != null) {
-          this.p = new double[original.p.GetLength(0), original.p.GetLength(1)];
-          Array.Copy(original.p, this.p, this.p.Length);
-        }
-        this.newData = new double[original.newData.GetLength(0), original.newData.GetLength(1)];
-        Array.Copy(original.newData, this.newData, this.newData.Length);
-        this.iter = original.iter;
-        this.currentMomentum = original.currentMomentum;
-        this.gains = new double[original.gains.GetLength(0), original.gains.GetLength(1)];
-        Array.Copy(original.gains, this.gains, this.gains.Length);
-        this.uY = new double[original.uY.GetLength(0), original.uY.GetLength(1)];
-        Array.Copy(original.uY, this.uY, this.uY.Length);
-        this.dY = new double[original.dY.GetLength(0), original.dY.GetLength(1)];
-        Array.Copy(original.dY, this.dY, this.dY.Length);
+        distance = cloner.Clone(original.distance);
+        random = cloner.Clone(original.random);
+        perplexity = original.perplexity;
+        exact = original.exact;
+        noDatapoints = original.noDatapoints;
+        finalMomentum = original.finalMomentum;
+        momSwitchIter = original.momSwitchIter;
+        stopLyingIter = original.stopLyingIter;
+        theta = original.theta;
+        eta = original.eta;
+        newDimensions = original.newDimensions;
+        if (original.valP != null) {
+          valP = new double[original.valP.Length];
+          Array.Copy(original.valP, valP, valP.Length);
+        }
+        if (original.rowP != null) {
+          rowP = new int[original.rowP.Length];
+          Array.Copy(original.rowP, rowP, rowP.Length);
+        }
+        if (original.colP != null) {
+          colP = new int[original.colP.Length];
+          Array.Copy(original.colP, colP, colP.Length);
+        }
+        if (original.p != null) {
+          p = new double[original.p.GetLength(0), original.p.GetLength(1)];
+          Array.Copy(original.p, p, p.Length);
+        }
+        newData = new double[original.newData.GetLength(0), original.newData.GetLength(1)];
+        Array.Copy(original.newData, newData, newData.Length);
+        iter = original.iter;
+        currentMomentum = original.currentMomentum;
+        gains = new double[original.gains.GetLength(0), original.gains.GetLength(1)];
+        Array.Copy(original.gains, gains, gains.Length);
+        uY = new double[original.uY.GetLength(0), original.uY.GetLength(1)];
+        Array.Copy(original.uY, uY, uY.Length);
+        dY = new double[original.dY.GetLength(0), original.dY.GetLength(1)];
+        Array.Copy(original.dY, dY, dY.Length);
       }
 
@@ -177 +178 @@
         this.stopLyingIter = stopLyingIter;
         this.momSwitchIter = momSwitchIter;
-        this.currentMomentum = momentum;
+        currentMomentum = momentum;
         this.finalMomentum = finalMomentum;
         this.eta = eta;
 
-
         // initialize
         noDatapoints = data.Length;
-        if(noDatapoints - 1 < 3 * perplexity)
+        if (noDatapoints - 1 < 3 * perplexity)
           throw new ArgumentException("Perplexity too large for the number of data points!");
 
@@ -192 +192 @@
         uY = new double[noDatapoints, newDimensions];
         gains = new double[noDatapoints, newDimensions];
-        for(var i = 0; i < noDatapoints; i++)
-          for(var j = 0; j < newDimensions; j++)
+        for (var i = 0; i < noDatapoints; i++)
+          for (var j = 0; j < newDimensions; j++)
             gains[i, j] = 1.0;
 
@@ -206 +206 @@
 
         // Lie about the P-values (factor is 4 in the MATLAB implementation)
-        if(exact) for(var i = 0; i < noDatapoints; i++) for(var j = 0; j < noDatapoints; j++) p[i, j] *= 12.0;
-        else for(var i = 0; i < rowP[noDatapoints]; i++) valP[i] *= 12.0;
+        if (exact) for (var i = 0; i < noDatapoints; i++) for (var j = 0; j < noDatapoints; j++) p[i, j] *= 12.0;
+        else for (var i = 0; i < rowP[noDatapoints]; i++) valP[i] *= 12.0;
 
         // Initialize solution (randomly)
         var rand = new NormalDistributedRandom(random, 0, 1);
-        for(var i = 0; i < noDatapoints; i++)
-          for(var j = 0; j < newDimensions; j++)
+        for (var i = 0; i < noDatapoints; i++)
+          for (var j = 0; j < newDimensions; j++)
             newData[i, j] = rand.NextDouble() * .0001;
       }
+      #endregion
 
       public double EvaluateError() {
@@ -222 +223 @@
       }
 
+      #region Helpers
       private static void CalculateApproximateSimilarities(T[] data, IDistance<T> distance, double perplexity, out int[] rowP, out int[] colP, out double[] valP) {
         // Compute asymmetric pairwise input similarities
@@ -233 +235 @@
         valP = sValP;
         var sumP = .0;
-        for(var i = 0; i < rowP[data.Length]; i++) sumP += valP[i];
-        for(var i = 0; i < rowP[data.Length]; i++) valP[i] /= sumP;
+        for (var i = 0; i < rowP[data.Length]; i++) sumP += valP[i];
+        for (var i = 0; i < rowP[data.Length]; i++) valP[i] /= sumP;
       }
 
@@ -242 +244 @@
         ComputeGaussianPerplexity(data, distance, p, perplexity);
         // Symmetrize input similarities
-        for(var n = 0; n < data.Length; n++) {
-          for(var m = n + 1; m < data.Length; m++) {
+        for (var n = 0; n < data.Length; n++) {
+          for (var m = n + 1; m < data.Length; m++) {
             p[n, m] += p[m, n];
             p[m, n] = p[n, m];
@@ -249 +251 @@
         }
         var sumP = .0;
-        for(var i = 0; i < data.Length; i++) for(var j = 0; j < data.Length; j++) sumP += p[i, j];
-        for(var i = 0; i < data.Length; i++) for(var j = 0; j < data.Length; j++) p[i, j] /= sumP;
+        for (var i = 0; i < data.Length; i++) for (var j = 0; j < data.Length; j++) sumP += p[i, j];
+        for (var i = 0; i < data.Length; i++) for (var j = 0; j < data.Length; j++) p[i, j] /= sumP;
         return p;
       }
 
       private static void ComputeGaussianPerplexity(IReadOnlyList<T> x, IDistance<T> distance, out int[] rowP, out int[] colP, out double[] valP, double perplexity, int k) {
-        if(perplexity > k) throw new ArgumentException("Perplexity should be lower than k!");
-
-        int n = x.Count;
+        if (perplexity > k) throw new ArgumentException("Perplexity should be lower than k!");
+
+        var n = x.Count;
         // Allocate the memory we need
         rowP = new int[n + 1];
@@ -264 +266 @@
         var curP = new double[n - 1];
         rowP[0] = 0;
-        for(var i = 0; i < n; i++) rowP[i + 1] = rowP[i] + k;
+        for (var i = 0; i < n; i++) rowP[i + 1] = rowP[i] + k;
 
         var objX = new List<IndexedItem<T>>();
-        for(var i = 0; i < n; i++) objX.Add(new IndexedItem<T>(i, x[i]));
+        for (var i = 0; i < n; i++) objX.Add(new IndexedItem<T>(i, x[i]));
 
         // Build ball tree on data set
@@ -273 +275 @@
 
         // Loop over all points to find nearest neighbors
-        for(var i = 0; i < n; i++) {
+        for (var i = 0; i < n; i++) {
           IList<IndexedItem<T>> indices;
           IList<double> distances;
@@ -285 +287 @@
           var minBeta = double.MinValue;
           var maxBeta = double.MaxValue;
-          const double tol = 1e-5;  
+          const double tol = 1e-5;
 
           // Iterate until we found a good perplexity
           var iter = 0; double sumP = 0;
-          while(!found && iter < 200) {  
+          while (!found && iter < 200) {
 
             // Compute Gaussian kernel row
-            for(var m = 0; m < k; m++) curP[m] = Math.Exp(-beta * distances[m + 1]); 
+            for (var m = 0; m < k; m++) curP[m] = Math.Exp(-beta * distances[m + 1]);
 
             // Compute entropy of current row
             sumP = double.Epsilon;
-            for(var m = 0; m < k; m++) sumP += curP[m];
+            for (var m = 0; m < k; m++) sumP += curP[m];
             var h = .0;
-            for(var m = 0; m < k; m++) h += beta * (distances[m + 1] * curP[m]); 
+            for (var m = 0; m < k; m++) h += beta * (distances[m + 1] * curP[m]);
             h = h / sumP + Math.Log(sumP);
 
             // Evaluate whether the entropy is within the tolerance level
             var hdiff = h - Math.Log(perplexity);
-            if(hdiff < tol && -hdiff < tol) {
+            if (hdiff < tol && -hdiff < tol) {
               found = true;
             } else {
-              if(hdiff > 0) {
+              if (hdiff > 0) {
                 minBeta = beta;
-                if(maxBeta.IsAlmost(double.MaxValue) || maxBeta.IsAlmost(double.MinValue))
+                if (maxBeta.IsAlmost(double.MaxValue) || maxBeta.IsAlmost(double.MinValue))
                   beta *= 2.0;
                 else
@@ -314 +316 @@
               } else {
                 maxBeta = beta;
-                if(minBeta.IsAlmost(double.MinValue) || minBeta.IsAlmost(double.MaxValue))
+                if (minBeta.IsAlmost(double.MinValue) || minBeta.IsAlmost(double.MaxValue))
                   beta /= 2.0;
                 else
@@ -326 +328 @@
 
           // Row-normalize current row of P and store in matrix
-          for(var m = 0; m < k; m++) curP[m] /= sumP;
-          for(var m = 0; m < k; m++) {
+          for (var m = 0; m < k; m++) curP[m] /= sumP;
+          for (var m = 0; m < k; m++) {
             colP[rowP[i] + m] = indices[m + 1].Index;
             valP[rowP[i] + m] = curP[m];
@@ -337 +339 @@
         var dd = ComputeDistances(x, distance);
 
-        int n = x.Length;
+        var n = x.Length;
         // Compute the Gaussian kernel row by row
-        for(var i = 0; i < n; i++) {
+        for (var i = 0; i < n; i++) {
           // Initialize some variables
           var found = false;
@@ -350 +352 @@
           // Iterate until we found a good perplexity
           var iter = 0;
-          while(!found && iter < 200) {      // 200 iterations as in tSNE implementation by van der Maarten
+          while (!found && iter < 200) {      // 200 iterations as in tSNE implementation by van der Maarten
 
             // Compute Gaussian kernel row
-            for(var m = 0; m < n; m++) p[i, m] = Math.Exp(-beta * dd[i][m]);
+            for (var m = 0; m < n; m++) p[i, m] = Math.Exp(-beta * dd[i][m]);
             p[i, i] = double.Epsilon;
 
             // Compute entropy of current row
             sumP = double.Epsilon;
-            for(var m = 0; m < n; m++) sumP += p[i, m];
+            for (var m = 0; m < n; m++) sumP += p[i, m];
             var h = 0.0;
-            for(var m = 0; m < n; m++) h += beta * (dd[i][m] * p[i, m]);
+            for (var m = 0; m < n; m++) h += beta * (dd[i][m] * p[i, m]);
             h = h / sumP + Math.Log(sumP);
 
             // Evaluate whether the entropy is within the tolerance level
             var hdiff = h - Math.Log(perplexity);
-            if(hdiff < tol && -hdiff < tol) {
+            if (hdiff < tol && -hdiff < tol) {
               found = true;
             } else {
-              if(hdiff > 0) {
+              if (hdiff > 0) {
                 minBeta = beta;
-                if(maxBeta.IsAlmost(double.MaxValue) || maxBeta.IsAlmost(double.MinValue))
+                if (maxBeta.IsAlmost(double.MaxValue) || maxBeta.IsAlmost(double.MinValue))
                   beta *= 2.0;
                 else
@@ -376 +378 @@
               } else {
                 maxBeta = beta;
-                if(minBeta.IsAlmost(double.MinValue) || minBeta.IsAlmost(double.MaxValue))
+                if (minBeta.IsAlmost(double.MinValue) || minBeta.IsAlmost(double.MaxValue))
                   beta /= 2.0;
                 else
@@ -388 +390 @@
 
           // Row normalize P
-          for(var m = 0; m < n; m++) p[i, m] /= sumP;
+          for (var m = 0; m < n; m++) p[i, m] /= sumP;
         }
       }
@@ -394 +396 @@
       private static double[][] ComputeDistances(T[] x, IDistance<T> distance) {
         var res = new double[x.Length][];
-        for(int r = 0; r < x.Length; r++) {
+        for (var r = 0; r < x.Length; r++) {
           var rowV = new double[x.Length];
           // all distances must be symmetric 
-          for(int c = 0; c < r; c++) {
+          for (var c = 0; c < r; c++) {
             rowV[c] = res[c][r];
           }
           rowV[r] = 0.0; // distance to self is zero for all distances
-          for(int c = r + 1; c < x.Length; c++) {
+          for (var c = r + 1; c < x.Length; c++) {
             rowV[c] = distance.Get(x[r], x[c]);
           }
@@ -418 +420 @@
         // Compute Q-matrix and normalization sum
         var sumQ = double.Epsilon;
-        for(var n1 = 0; n1 < n; n1++) {
-          for(var m = 0; m < n; m++) {
-            if(n1 != m) {
+        for (var n1 = 0; n1 < n; n1++) {
+          for (var m = 0; m < n; m++) {
+            if (n1 != m) {
               q[n1, m] = 1 / (1 + dd[n1, m]);
               sumQ += q[n1, m];
@@ -426 +428 @@
           }
         }
-        for(var i = 0; i < n; i++) for(var j = 0; j < n; j++) q[i, j] /= sumQ;
+        for (var i = 0; i < n; i++) for (var j = 0; j < n; j++) q[i, j] /= sumQ;
 
         // Sum t-SNE error
         var c = .0;
-        for(var i = 0; i < n; i++)
-          for(var j = 0; j < n; j++) {
+        for (var i = 0; i < n; i++)
+          for (var j = 0; j < n; j++) {
             c += p[i, j] * Math.Log((p[i, j] + float.Epsilon) / (q[i, j] + float.Epsilon));
           }
@@ -437 +439 @@
       }
 
-      // The mapping of the approximate tSNE looks good but the error curve never changes.
       private static double EvaluateErrorApproximate(IReadOnlyList<int> rowP, IReadOnlyList<int> colP, IReadOnlyList<double> valP, double[,] y, double theta) {
         // Get estimate of normalization term
@@ -444 +445 @@
         var tree = new SpacePartitioningTree(y);
         var buff = new double[d];
-        double sumQ = 0.0;
-        for(var i = 0; i < n; i++) tree.ComputeNonEdgeForces(i, theta, buff, ref sumQ);
+        var sumQ = 0.0;
+        for (var i = 0; i < n; i++) tree.ComputeNonEdgeForces(i, theta, buff, ref sumQ);
 
         // Loop over all edges to compute t-SNE error
         var c = .0;
-        for(var k = 0; k < n; k++) {
-          for(var i = rowP[k]; i < rowP[k + 1]; i++) {
+        for (var k = 0; k < n; k++) {
+          for (var i = rowP[k]; i < rowP[k + 1]; i++) {
             var q = .0;
-            for(var j = 0; j < d; j++) buff[j] = y[k, j];
-            for(var j = 0; j < d; j++) buff[j] -= y[colP[i], j];
-            for(var j = 0; j < d; j++) q += buff[j] * buff[j];
+            for (var j = 0; j < d; j++) buff[j] = y[k, j];
+            for (var j = 0; j < d; j++) buff[j] -= y[colP[i], j];
+            for (var j = 0; j < d; j++) q += buff[j] * buff[j];
             q = (1.0 / (1.0 + q)) / sumQ;
             c += valP[i] * Math.Log((valP[i] + float.Epsilon) / (q + float.Epsilon));
@@ -466 +467 @@
         var n = rowP.Count - 1;
         var rowCounts = new int[n];
-        for(var j = 0; j < n; j++) {
-          for(var i = rowP[j]; i < rowP[j + 1]; i++) {
+        for (var j = 0; j < n; j++) {
+          for (var i = rowP[j]; i < rowP[j + 1]; i++) {
 
             // Check whether element (col_P[i], n) is present
             var present = false;
-            for(var m = rowP[colP[i]]; m < rowP[colP[i] + 1]; m++) {
-              if(colP[m] == j) present = true;
+            for (var m = rowP[colP[i]]; m < rowP[colP[i] + 1]; m++) {
+              if (colP[m] == j) present = true;
             }
-            if(present) rowCounts[j]++;
+            if (present) rowCounts[j]++;
             else {
               rowCounts[j]++;
@@ -482 +483 @@
         }
         var noElem = 0;
-        for(var i = 0; i < n; i++) noElem += rowCounts[i];
+        for (var i = 0; i < n; i++) noElem += rowCounts[i];
 
         // Allocate memory for symmetrized matrix
@@ -491 +492 @@
         // Construct new row indices for symmetric matrix
         symRowP[0] = 0;
-        for(var i = 0; i < n; i++) symRowP[i + 1] = symRowP[i] + rowCounts[i];
+        for (var i = 0; i < n; i++) symRowP[i + 1] = symRowP[i] + rowCounts[i];
 
         // Fill the result matrix
         var offset = new int[n];
-        for(var j = 0; j < n; j++) {
-          for(var i = rowP[j]; i < rowP[j + 1]; i++) {                                  // considering element(n, colP[i])
+        for (var j = 0; j < n; j++) {
+          for (var i = rowP[j]; i < rowP[j + 1]; i++) {                                  // considering element(n, colP[i])
 
             // Check whether element (col_P[i], n) is present
             var present = false;
-            for(var m = rowP[colP[i]]; m < rowP[colP[i] + 1]; m++) {
-              if(colP[m] != j) continue;
+            for (var m = rowP[colP[i]]; m < rowP[colP[i] + 1]; m++) {
+              if (colP[m] != j) continue;
               present = true;
-              if(j > colP[i]) continue; // make sure we do not add elements twice
+              if (j > colP[i]) continue; // make sure we do not add elements twice
               symColP[symRowP[j] + offset[j]] = colP[i];
               symColP[symRowP[colP[i]] + offset[colP[i]]] = j;
@@ -511 +512 @@
 
             // If (colP[i], n) is not present, there is no addition involved
-            if(!present) {
+            if (!present) {
               symColP[symRowP[j] + offset[j]] = colP[i];
               symColP[symRowP[colP[i]] + offset[colP[i]]] = j;
@@ -519 +520 @@
 
             // Update offsets
-            if(present && (j > colP[i])) continue;
+            if (present && (j > colP[i])) continue;
             offset[j]++;
-            if(colP[i] != j) offset[colP[i]]++;
-          }
-        }
-
-        for(var i = 0; i < noElem; i++) symValP[i] /= 2.0;
-      }
+            if (colP[i] != j) offset[colP[i]]++;
+          }
+        }
+
+        for (var i = 0; i < noElem; i++) symValP[i] /= 2.0;
+      }
+      #endregion
     }
 
     /// <summary>
-    /// Simple interface to tSNE
+    /// Static interface to tSNE
     /// </summary>
     /// <param name="data"></param>
@@ -548 +550 @@
       int newDimensions = 2, double perplexity = 25, int iterations = 1000,
       double theta = 0,
-      int stopLyingIter = 250, int momSwitchIter = 250, double momentum = .5,
-      double finalMomentum = .8, double eta = 200.0
+      int stopLyingIter = 0, int momSwitchIter = 0, double momentum = .5,
+      double finalMomentum = .8, double eta = 10.0
       ) {
       var state = CreateState(data, distance, random, newDimensions, perplexity,
         theta, stopLyingIter, momSwitchIter, momentum, finalMomentum, eta);
 
-      for(int i = 0; i < iterations - 1; i++) {
+      for (var i = 0; i < iterations - 1; i++) {
         Iterate(state);
       }
@@ -562 +564 @@
     public static TSNEState CreateState(T[] data, IDistance<T> distance, IRandom random,
       int newDimensions = 2, double perplexity = 25, double theta = 0,
-      int stopLyingIter = 250, int momSwitchIter = 250, double momentum = .5,
-      double finalMomentum = .8, double eta = 200.0
+      int stopLyingIter = 0, int momSwitchIter = 0, double momentum = .5,
+      double finalMomentum = .8, double eta = 10.0
       ) {
       return new TSNEState(data, distance, random, newDimensions, perplexity, theta, stopLyingIter, momSwitchIter, momentum, finalMomentum, eta);
     }
 
-
     public static double[,] Iterate(TSNEState state) {
-      if(state.exact)
+      if (state.exact)
         ComputeExactGradient(state.p, state.newData, state.noDatapoints, state.newDimensions, state.dY);
       else
@@ -576 +577 @@
 
       // Update gains
-      for(var i = 0; i < state.noDatapoints; i++) {
-        for(var j = 0; j < state.newDimensions; j++) {
+      for (var i = 0; i < state.noDatapoints; i++) {
+        for (var j = 0; j < state.newDimensions; j++) {
           state.gains[i, j] = Math.Sign(state.dY[i, j]) != Math.Sign(state.uY[i, j])
             ? state.gains[i, j] + .2  // +0.2 nd *0.8 are used in two separate implementations of tSNE -> seems to be correct
             : state.gains[i, j] * .8;
 
-          if(state.gains[i, j] < .01) state.gains[i, j] = .01;
+          if (state.gains[i, j] < .01) state.gains[i, j] = .01;
         }
       }
@@ -588 +589 @@
 
       // Perform gradient update (with momentum and gains)
-      for(var i = 0; i < state.noDatapoints; i++)
-        for(var j = 0; j < state.newDimensions; j++)
+      for (var i = 0; i < state.noDatapoints; i++)
+        for (var j = 0; j < state.newDimensions; j++)
           state.uY[i, j] = state.currentMomentum * state.uY[i, j] - state.eta * state.gains[i, j] * state.dY[i, j];
 
-      for(var i = 0; i < state.noDatapoints; i++)
-        for(var j = 0; j < state.newDimensions; j++)
+      for (var i = 0; i < state.noDatapoints; i++)
+        for (var j = 0; j < state.newDimensions; j++)
           state.newData[i, j] = state.newData[i, j] + state.uY[i, j];
 
@@ -600 +601 @@
 
       // Stop lying about the P-values after a while, and switch momentum
-      if(state.iter == state.stopLyingIter) {
-        if(state.exact)
-          for(var i = 0; i < state.noDatapoints; i++)
-            for(var j = 0; j < state.noDatapoints; j++)
+      if (state.iter == state.stopLyingIter) {
+        if (state.exact)
+          for (var i = 0; i < state.noDatapoints; i++)
+            for (var j = 0; j < state.noDatapoints; j++)
               state.p[i, j] /= 12.0;
         else
-          for(var i = 0; i < state.rowP[state.noDatapoints]; i++)
+          for (var i = 0; i < state.rowP[state.noDatapoints]; i++)
             state.valP[i] /= 12.0;
       }
 
-      if(state.iter == state.momSwitchIter)
+      if (state.iter == state.momSwitchIter)
         state.currentMomentum = state.finalMomentum;
 
@@ -617 +618 @@
     }
 
-
-
+    #region Helpers
     private static void ComputeApproximateGradient(int[] rowP, int[] colP, double[] valP, double[,] y, int n, int d, double[,] dC, double theta) {
       var tree = new SpacePartitioningTree(y);
-      double sumQ = 0.0;
+      var sumQ = 0.0;
       var posF = new double[n, d];
       var negF = new double[n, d];
-      tree.ComputeEdgeForces(rowP, colP, valP, n, posF);
+      SpacePartitioningTree.ComputeEdgeForces(rowP, colP, valP, n, posF, y, d);
       var row = new double[d];
-      for(var n1 = 0; n1 < n; n1++) {
-        Array.Clear(row,0,row.Length);
+      for (var n1 = 0; n1 < n; n1++) {
+        Array.Clear(row, 0, row.Length);
         tree.ComputeNonEdgeForces(n1, theta, row, ref sumQ);
-        Buffer.BlockCopy(row, 0, negF, (sizeof(double) * n1 * d), d*sizeof(double));
+        Buffer.BlockCopy(row, 0, negF, (sizeof(double) * n1 * d), d * sizeof(double));
       }
 
       // Compute final t-SNE gradient
       for (var i = 0; i < n; i++)
-        for(var j = 0; j < d; j++) {
+        for (var j = 0; j < d; j++) {
           dC[i, j] = posF[i, j] - negF[i, j] / sumQ;
         }
@@ -640 +640 @@
 
     private static void ComputeExactGradient(double[,] p, double[,] y, int n, int d, double[,] dC) {
-
       // Make sure the current gradient contains zeros
-      for(var i = 0; i < n; i++) for(var j = 0; j < d; j++) dC[i, j] = 0.0;
+      for (var i = 0; i < n; i++) for (var j = 0; j < d; j++) dC[i, j] = 0.0;
 
       // Compute the squared Euclidean distance matrix
@@ -651 +650 @@
       var q = new double[n, n];
       var sumQ = .0;
-      for(var n1 = 0; n1 < n; n1++) {
-        for(var m = 0; m < n; m++) {
-          if(n1 == m) continue;
+      for (var n1 = 0; n1 < n; n1++) {
+        for (var m = 0; m < n; m++) {
+          if (n1 == m) continue;
           q[n1, m] = 1 / (1 + dd[n1, m]);
           sumQ += q[n1, m];
@@ -660 +659 @@
 
       // Perform the computation of the gradient
-      for(var n1 = 0; n1 < n; n1++) {
-        for(var m = 0; m < n; m++) {
-          if(n1 == m) continue;
+      for (var n1 = 0; n1 < n; n1++) {
+        for (var m = 0; m < n; m++) {
+          if (n1 == m) continue;
           var mult = (p[n1, m] - q[n1, m] / sumQ) * q[n1, m];
-          for(var d1 = 0; d1 < d; d1++) {
+          for (var d1 = 0; d1 < d; d1++) {
             dC[n1, d1] += (y[n1, d1] - y[m, d1]) * mult;
           }
@@ -673 +672 @@
     private static void ComputeSquaredEuclideanDistance(double[,] x, int n, int d, double[,] dd) {
       var dataSums = new double[n];
-      for(var i = 0; i < n; i++) {
-        for(var j = 0; j < d; j++) {
+      for (var i = 0; i < n; i++) {
+        for (var j = 0; j < d; j++) {
           dataSums[i] += x[i, j] * x[i, j];
         }
       }
-      for(var i = 0; i < n; i++) {
-        for(var m = 0; m < n; m++) {
+      for (var i = 0; i < n; i++) {
+        for (var m = 0; m < n; m++) {
           dd[i, m] = dataSums[i] + dataSums[m];
         }
       }
-      for(var i = 0; i < n; i++) {
+      for (var i = 0; i < n; i++) {
         dd[i, i] = 0.0;
-        for(var m = i + 1; m < n; m++) {
+        for (var m = i + 1; m < n; m++) {
           dd[i, m] = 0.0;
-          for(var j = 0; j < d; j++) {
+          for (var j = 0; j < d; j++) {
             dd[i, m] += (x[i, j] - x[m, j]) * (x[i, j] - x[m, j]);
           }
@@ -700 +699 @@
       var d = x.GetLength(1);
       var mean = new double[d];
-      for(var i = 0; i < n; i++) {
-        for(var j = 0; j < d; j++) {
+      for (var i = 0; i < n; i++) {
+        for (var j = 0; j < d; j++) {
           mean[j] += x[i, j];
         }
       }
-      for(var i = 0; i < d; i++) {
+      for (var i = 0; i < d; i++) {
         mean[i] /= n;
       }
       // Subtract data mean
-      for(var i = 0; i < n; i++) {
-        for(var j = 0; j < d; j++) {
+      for (var i = 0; i < n; i++) {
+        for (var j = 0; j < d; j++) {
           x[i, j] -= mean[j];
         }
       }
     }
+    #endregion
   }
 }

trunk/sources/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/VantagePointTree.cs

@@ -68 +68 @@
   /// </summary>
   /// <typeparam name="T"></typeparam>
-  public class VantagePointTree<T> : IVantagePointTree<T> {
-    #region properties
+  public class VantagePointTree<T> {
     private readonly List<T> items;
     private readonly Node root;
     private readonly IDistance<T> distance;
-    #endregion
 
     public VantagePointTree(IDistance<T> distance, IEnumerable<T> items) : base() {
@@ -82 +80 @@
     }
 
+    /// <summary>
+    /// provides the k-nearest neighbours to a certain target element
+    /// </summary>
+    /// <param name="target">The target element</param>
+    /// <param name="k">How many neighbours</param>
+    /// <param name="results">The nearest neighbouring elements</param>
+    /// <param name="distances">The distances form the target corresponding to the neighbouring elements</param>
     public void Search(T target, int k, out IList<T> results, out IList<double> distances) {
       var heap = new PriorityQueue<double, IndexedItem<double>>(double.MaxValue, double.MinValue, k);
-      double tau = double.MaxValue;
+      var tau = double.MaxValue;
       Search(root, target, k, heap, ref tau);
       var res = new List<T>();
       var dist = new List<double>();
       while (heap.Size > 0) {
-        res.Add(items[heap.PeekMinValue().Index]);          // actually max distance
+        res.Add(items[heap.PeekMinValue().Index]);// actually max distance
         dist.Add(heap.PeekMinValue().Value);
         heap.RemoveMin();
       }
-      res.Reverse();  
+      res.Reverse();
       dist.Reverse();
       results = res;
@@ -104 +109 @@
       if (dist < tau) {
         if (heap.Size == k) heap.RemoveMin();   // remove furthest node from result list (if we already have k results) 
-        heap.Insert(-dist, new IndexedItem<double>(node.index, dist));      
+        heap.Insert(-dist, new IndexedItem<double>(node.index, dist));
         if (heap.Size == k) tau = heap.PeekMinValue().Value;
       }