Index: branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/Interfaces/TSNEInterfaces/ISpacePartitioningTree.cs
===================================================================
--- branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/Interfaces/TSNEInterfaces/ISpacePartitioningTree.cs	(revision 14788)
+++ branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/Interfaces/TSNEInterfaces/ISpacePartitioningTree.cs	(revision 14806)
@@ -25,7 +25,5 @@
   public interface ISpacePartitioningTree {
 
-    double[,] Data { set; }
     int GetDepth();
-    ISpacePartitioningTree GetParent();
     bool Insert(int newIndex);
     void Subdivide();
Index: branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/SpacePartitioningTree.cs
===================================================================
--- branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/SpacePartitioningTree.cs	(revision 14788)
+++ branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/SpacePartitioningTree.cs	(revision 14806)
@@ -77,5 +77,5 @@
 
     // Indices in this space-partitioning tree node, corresponding center-of-mass, and list of all children
-    public double[,] Data { private get; set; } // TODO public setter, private getter?
+    private double[,] data;
 
     private double[] centerOfMass;
@@ -91,17 +91,17 @@
       var meanY = new double[d];
       var minY = new double[d];
-      for (var i = 0; i < d; i++) minY[i] = double.MaxValue;
+      for(var i = 0; i < d; i++) minY[i] = double.MaxValue;
       var maxY = new double[d];
-      for (var i = 0; i < d; i++) maxY[i] = double.MinValue;
-      for (uint i = 0; i < n; i++) {
-        for (uint j = 0; j < d; j++) {
+      for(var i = 0; i < d; i++) maxY[i] = double.MinValue;
+      for(uint i = 0; i < n; i++) {
+        for(uint j = 0; j < d; j++) {
           meanY[j] = inpData[i, j];
-          if (inpData[i, j] < minY[j]) minY[j] = inpData[i, j];
-          if (inpData[i, j] > maxY[j]) maxY[j] = inpData[i, j];
-        }
-      }
-      for (var i = 0; i < d; i++) meanY[i] /= n;
+          if(inpData[i, j] < minY[j]) minY[j] = inpData[i, j];
+          if(inpData[i, j] > maxY[j]) maxY[j] = inpData[i, j];
+        }
+      }
+      for(var i = 0; i < d; i++) meanY[i] /= n;
       var width = new double[d];
-      for (var i = 0; i < d; i++) width[i] = Math.Max(maxY[i] - meanY[i], meanY[i] - maxY[i]) + 1e-5;
+      for(var i = 0; i < d; i++) width[i] = Math.Max(maxY[i] - meanY[i], meanY[i] - maxY[i]) + 1e-5; // TODO constant?
       Init(null, inpData, meanY, width);
       Fill(n);
@@ -122,15 +122,15 @@
       // Ignore objects which do not belong in this quad tree
       var point = new double[dimension];
-      Buffer.BlockCopy(Data, sizeof(double) * dimension * newIndex, point, 0, sizeof(double) * dimension);
-      if (!boundary.ContainsPoint(point)) return false;
+      Buffer.BlockCopy(data, sizeof(double) * dimension * newIndex, point, 0, sizeof(double) * dimension);
+      if(!boundary.ContainsPoint(point)) return false;
       cumulativeSize++;
       // Online update of cumulative size and center-of-mass
       var mult1 = (double)(cumulativeSize - 1) / cumulativeSize;
       var mult2 = 1.0 / cumulativeSize;
-      for (var i = 0; i < dimension; i++) centerOfMass[i] *= mult1;
-      for (var i = 0; i < dimension; i++) centerOfMass[i] += mult2 * point[i];
+      for(var i = 0; i < dimension; i++) centerOfMass[i] *= mult1;
+      for(var i = 0; i < dimension; i++) centerOfMass[i] += mult2 * point[i];
 
       // 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 < QT_NODE_CAPACITY) {
         index[size] = newIndex;
         size++;
@@ -140,19 +140,19 @@
       // Don't add duplicates for now (this is not very nice)
       var anyDuplicate = false;
-      for (uint n = 0; n < size; n++) {
+      for(uint n = 0; n < size; n++) {
         var duplicate = true;
-        for (var d = 0; d < dimension; d++) {
-          if (Math.Abs(point[d] - Data[index[n], d]) < double.Epsilon) continue;
+        for(var d = 0; d < dimension; d++) {
+          if(Math.Abs(point[d] - data[index[n], d]) < double.Epsilon) continue;
           duplicate = false; break;
         }
         anyDuplicate = anyDuplicate | duplicate;
       }
-      if (anyDuplicate) return true;
+      if(anyDuplicate) return true;
 
       // Otherwise, we need to subdivide the current cell
-      if (isLeaf) Subdivide();
+      if(isLeaf) Subdivide();
       // Find out where the point can be inserted
-      for (var i = 0; i < noChildren; i++) {
-        if (children[i].Insert(newIndex)) return true;
+      for(var i = 0; i < noChildren; i++) {
+        if(children[i].Insert(newIndex)) return true;
       }
 
@@ -165,20 +165,20 @@
       var newCorner = new double[dimension];
       var newWidth = new double[dimension];
-      for (var i = 0; i < noChildren; i++) {
+      for(var i = 0; i < noChildren; i++) {
         var div = 1;
-        for (var d = 0; d < dimension; d++) {
+        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);
+          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);
+        children[i] = new SpacePartitioningTree(this, data, newCorner, newWidth);
       }
 
       // Move existing points to correct children
-      for (var i = 0; i < size; i++) {
+      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]);
+        for(var j = 0; j < noChildren; j++) {
+          if(!success) success = children[j].Insert(index[i]);
         }
         index[i] = int.MaxValue;
@@ -192,9 +192,9 @@
     public bool IsCorrect() {
       var row = new double[dimension];
-      for (var n = 0; n < size; n++) Buffer.BlockCopy(Data, sizeof(double) * dimension * n, row, 0, sizeof(double) * dimension);
-      if (!boundary.ContainsPoint(row)) return false;
-      if (isLeaf) return true;
+      for(var n = 0; n < size; n++) Buffer.BlockCopy(data, sizeof(double) * dimension * 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();
+      for(var i = 0; i < noChildren; i++) correct = correct && children[i].IsCorrect();
       return correct;
     }
@@ -206,9 +206,9 @@
     public int GetAllIndices(int[] indices, int loc) {
       // Gather indices in current quadrant
-      for (var i = 0; i < size; i++) indices[loc + i] = index[i];
+      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);
+      if(isLeaf) return loc;
+      for(var i = 0; i < noChildren; i++) loc = children[i].GetAllIndices(indices, loc);
       return loc;
     }
@@ -218,21 +218,21 @@
     }
 
-    public void ComputeNonEdgeForces(int pointIndex, double theta, double[] negF, ref double sumQ)
-    {
+    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
-      if (cumulativeSize == 0 || (isLeaf && size == 1 && index[0] == pointIndex)) return;
+      if(cumulativeSize == 0 || (isLeaf && size == 1 && index[0] == pointIndex)) return;
 
       // Compute distance between point and center-of-mass
+      // TODO: squared distance with normalized axes is used here!
       var D = .0;
-      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];
+      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];
 
       // Check whether we can use this node as a "summary"
       var maxWidth = 0.0;
-      for (var d = 0; d < dimension; d++) {
+      for(var d = 0; d < dimension; d++) {
         var curWidth = boundary.GetWidth(d);
         maxWidth = (maxWidth > curWidth) ? maxWidth : curWidth;
       }
-      if (isLeaf || maxWidth / Math.Sqrt(D) < theta) {
+      if(isLeaf || maxWidth / Math.Sqrt(D) < theta) {
 
         // Compute and add t-SNE force between point and current node
@@ -241,9 +241,9 @@
         sumQ += mult;
         mult *= D;
-        for (var d = 0; d < dimension; d++) negF[d] += mult * buff[d];
+        for(var d = 0; d < dimension; d++) negF[d] += mult * buff[d];
       } else {
 
         // Recursively apply Barnes-Hut to children
-        for (var i = 0; i < noChildren; i++) children[i].ComputeNonEdgeForces(pointIndex, theta, negF, ref sumQ);
+        for(var i = 0; i < noChildren; i++) children[i].ComputeNonEdgeForces(pointIndex, theta, negF, ref sumQ);
       }
     }
@@ -251,15 +251,16 @@
     public void ComputeEdgeForces(int[] rowP, int[] colP, double[] valP, int n, double[,] posF) {
       // Loop over all edges in the graph
-      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++) {
 
           // Compute pairwise distance and Q-value
+          // uses squared distance
           var d = 1.0;
-          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];
+          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];
           d = valP[i] / d;
 
           // Sum positive force
-          for (var j = 0; j < dimension; j++) posF[k, j] += d * buff[j];
+          for(var j = 0; j < dimension; j++) posF[k, j] += d * buff[j];
         }
       }
@@ -268,5 +269,5 @@
     #region Helpers
     private void Fill(int n) {
-      for (var i = 0; i < n; i++) Insert(i);
+      for(var i = 0; i < n; i++) Insert(i);
     }
     private void Init(SpacePartitioningTree p, double[,] inpData, IEnumerable<double> inpCorner, IEnumerable<double> inpWidth) {
@@ -274,6 +275,6 @@
       dimension = inpData.GetLength(1);
       noChildren = 2;
-      for (uint i = 1; i < dimension; i++) noChildren *= 2;
-      Data = inpData;
+      for(uint i = 1; i < dimension; i++) noChildren *= 2;
+      data = inpData;
       isLeaf = true;
       size = 0;
@@ -291,5 +292,5 @@
 
 
-    private class Cell  {
+    private class Cell {
       private readonly uint dimension;
       private readonly double[] corner;
@@ -315,6 +316,6 @@
       }
       public bool ContainsPoint(double[] point) {
-        for (var d = 0; d < dimension; d++)
-          if (corner[d] - width[d] > point[d] || corner[d] + width[d] < point[d]) return false;
+        for(var d = 0; d < dimension; d++)
+          if(corner[d] - width[d] > point[d] || corner[d] + width[d] < point[d]) return false;
         return true;
       }
Index: branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/TSNEStatic.cs
===================================================================
--- branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/TSNEStatic.cs	(revision 14788)
+++ branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/TSNEStatic.cs	(revision 14806)
@@ -70,37 +70,96 @@
     public sealed class TSNEState : DeepCloneable {
       // initialized once
+      [Storable]
       public IDistance<T> distance;
+      [Storable]
       public IRandom random;
+      [Storable]
       public double perplexity;
+      [Storable]
       public bool exact;
+      [Storable]
       public int noDatapoints;
+      [Storable]
       public double finalMomentum;
+      [Storable]
       public int momSwitchIter;
+      [Storable]
       public int stopLyingIter;
+      [Storable]
       public double theta;
+      [Storable]
       public double eta;
+      [Storable]
       public int newDimensions;
 
       // for approximate version: sparse representation of similarity/distance matrix
+      [Storable]
       public double[] valP; // similarity/distance
+      [Storable]
       public int[] rowP; // row index
+      [Storable]
       public int[] colP; // col index
 
       // for exact version: dense representation of distance/similarity matrix
+      [Storable]
       public double[,] p;
 
       // mapped data
+      [Storable]
       public double[,] newData;
 
+      [Storable]
       public int iter;
+      [Storable]
       public double currentMomentum;
 
       // helper variables (updated in each iteration)
+      [Storable]
       public double[,] gains;
+      [Storable]
       public double[,] uY;
+      [Storable]
       public double[,] dY;
 
       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);
+      }
+
       public override IDeepCloneable Clone(Cloner cloner) {
         return new TSNEState(this, cloner);
@@ -122,5 +181,6 @@
         // initialize
         noDatapoints = data.Length;
-        if (noDatapoints - 1 < 3 * perplexity) throw new ArgumentException("Perplexity too large for the number of data points!");
+        if(noDatapoints - 1 < 3 * perplexity)
+          throw new ArgumentException("Perplexity too large for the number of data points!");
 
         exact = Math.Abs(theta) < double.Epsilon;
@@ -129,6 +189,6 @@
         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;
 
@@ -139,20 +199,22 @@
 
         //Calculate Similarities
-        if (exact) p = CalculateExactSimilarites(data, distance, perplexity);
+        if(exact) p = CalculateExactSimilarites(data, distance, perplexity);
         else CalculateApproximateSimilarities(data, distance, perplexity, out rowP, out colP, out valP);
 
         // Lie about the P-values
-        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++)
-            newData[i, j] = rand.NextDouble() * .0001;  // TODO const
+        for(var i = 0; i < noDatapoints; i++)
+          for(var j = 0; j < newDimensions; j++)
+            newData[i, j] = rand.NextDouble() * .0001;  // TODO const  ?
       }
 
       public double EvaluateError() {
-        return exact ? EvaluateErrorExact(p, newData, noDatapoints, newDimensions) : EvaluateErrorApproximate(rowP, colP, valP, newData, theta);
+        return exact ?
+          EvaluateErrorExact(p, newData, noDatapoints, newDimensions) :
+          EvaluateErrorApproximate(rowP, colP, valP, newData, theta);
       }
 
@@ -168,7 +230,8 @@
         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;
+      }
+
       private static double[,] CalculateExactSimilarites(T[] data, IDistance<T> distance, double perplexity) {
         // Compute similarities
@@ -176,6 +239,6 @@
         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];
@@ -183,11 +246,11 @@
         }
         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!");
+        if(perplexity > k) throw new ArgumentException("Perplexity should be lower than k!");
 
         int n = x.Count;
@@ -198,8 +261,8 @@
         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
@@ -207,5 +270,5 @@
 
         // 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;
@@ -223,24 +286,24 @@
           // Iterate until we found a good perplexity
           var iter = 0; double sumP = 0;
-          while (!found && iter < 200) {
+          while(!found && iter < 200) {  // TODO 200 iterations always ok?
 
             // 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]); // TODO 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]); // TODO: distances m+1?
             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
@@ -248,5 +311,5 @@
               } 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
@@ -260,6 +323,6 @@
 
           // 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];
@@ -273,5 +336,5 @@
         int 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;
@@ -284,25 +347,25 @@
           // Iterate until we found a good perplexity
           var iter = 0;
-          while (!found && iter < 200) {       // TODO constant
+          while(!found && iter < 200) {       // TODO constant
 
             // 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
@@ -310,5 +373,5 @@
               } 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
@@ -322,13 +385,25 @@
 
           // 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;
         }
       }
 
       private static double[][] ComputeDistances(T[] x, IDistance<T> distance) {
-        return x.Select(m => x.Select(n => distance.Get(m, n)).ToArray()).ToArray();
-      }
-
-
+        var res = new double[x.Length][];
+        for(int r = 0; r < x.Length; r++) {
+          var rowV = new double[x.Length];
+          // all distances must be symmetric 
+          for(int 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++) {
+            rowV[c] = distance.Get(x[r], x[c]);
+          }
+          res[r] = rowV;
+        }
+        return res;
+        // return x.Select(m => x.Select(n => distance.Get(m, n)).ToArray()).ToArray();
+      }
 
       private static double EvaluateErrorExact(double[,] p, double[,] y, int n, int d) {
@@ -336,11 +411,11 @@
         var dd = new double[n, n];
         var q = new double[n, n];
-        ComputeSquaredEuclideanDistance(y, n, d, dd);
+        ComputeSquaredEuclideanDistance(y, n, d, dd); // TODO: we use Euclidian distance regardless of the actual distance function
 
         // 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];
@@ -348,14 +423,17 @@
           }
         }
-        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));
           }
         return c;
       }
+
+      // TODO: there seems to be a bug in the error approximation.
+      // 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
@@ -365,14 +443,14 @@
         var buff = new double[d];
         double sumQ = 0.0;
-        for (var i = 0; i < n; i++) tree.ComputeNonEdgeForces(i, theta, buff, ref sumQ);
+        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];     // TODO: squared error is used here!
             q = 1.0 / (1.0 + q) / sumQ;
             c += valP[i] * Math.Log((valP[i] + float.Epsilon) / (q + float.Epsilon));
@@ -386,13 +464,13 @@
         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]++;
@@ -402,5 +480,5 @@
         }
         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
@@ -411,17 +489,17 @@
         // 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;
@@ -431,5 +509,5 @@
 
             // 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;
@@ -439,12 +517,11 @@
 
             // Update offsets
-            if (present && (j > colP[i])) continue;
+            if(present && (j > colP[i])) continue;
             offset[j]++;
-            if (colP[i] != j) offset[colP[i]]++;
-          }
-        }
-
-        // Divide the result by two
-        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;
       }
 
@@ -459,5 +536,5 @@
 
     public static double[,] Iterate(TSNEState state) {
-      if (state.exact)
+      if(state.exact)
         ComputeExactGradient(state.p, state.newData, state.noDatapoints, state.newDimensions, state.dY);
       else
@@ -465,11 +542,11 @@
 
       // 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
             : state.gains[i, j] * .8; // 20% up or 20% down // TODO: +0.2?!
 
-          if (state.gains[i, j] < .01) state.gains[i, j] = .01;
+          if(state.gains[i, j] < .01) state.gains[i, j] = .01; // TODO why limit the gains?
         }
       }
@@ -477,10 +554,10 @@
 
       // 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];
 
@@ -489,12 +566,12 @@
       // 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++) state.p[i, j] /= 12.0;                                   //XXX why 12?
+      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;                                   //XXX why 12?
         else
-          for (var i = 0; i < state.rowP[state.noDatapoints]; i++) state.valP[i] /= 12.0;                       // XXX are we not scaling all values? 
-      }
-
-      if (state.iter == state.momSwitchIter)
+          for(var i = 0; i < state.rowP[state.noDatapoints]; i++) state.valP[i] /= 12.0;                       // XXX are we not scaling all values? 
+      }
+
+      if(state.iter == state.momSwitchIter)
         state.currentMomentum = state.finalMomentum;
 
@@ -511,5 +588,5 @@
       tree.ComputeEdgeForces(rowP, colP, valP, n, posF);
       var row = new double[d];
-      for (var n1 = 0; n1 < n; n1++) {
+      for(var n1 = 0; n1 < n; n1++) {
         Buffer.BlockCopy(negF, (sizeof(double) * n1 * d), row, 0, d);
         tree.ComputeNonEdgeForces(n1, theta, row, ref sumQ);
@@ -517,6 +594,6 @@
 
       // Compute final t-SNE gradient
-      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++) {
           dC[i, j] = posF[i, j] - negF[i, j] / sumQ;
         }
@@ -526,16 +603,16 @@
 
       // 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
       var dd = new double[n, n];
-      ComputeSquaredEuclideanDistance(y, n, d, dd);
+      ComputeSquaredEuclideanDistance(y, n, d, dd); // TODO: we use Euclidian distance regardless which distance function is actually set!
 
       // Compute Q-matrix and normalization sum
       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];
@@ -544,9 +621,9 @@
 
       // 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;
           }
@@ -557,19 +634,19 @@
     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]);
           }
@@ -584,15 +661,15 @@
       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];
         }
