Changeset 14855 for branches

Timestamp:

04/12/17 12:15:37 (7 years ago)

Author:

gkronber

Message:

#2700: made some changes / bug-fixes while reviewing

Location:

branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE

Files:

• SpacePartitioningTree.cs (modified) (14 diffs)
• TSNEAlgorithm.cs (modified) (1 diff)
• TSNEStatic.cs (modified) (2 diffs)
• VantagePointTree.cs (modified) (4 diffs)

branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/SpacePartitioningTree.cs

@@ -76 +76 @@
     private Cell boundary;
 
+    private double[,] data;
+
     // Indices in this space-partitioning tree node, corresponding center-of-mass, and list of all children
-    private double[,] data;
-
     private double[] centerOfMass;
     private readonly int[] index = new int[QT_NODE_CAPACITY];
@@ -91 +91 @@
       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++) {
-          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;
+      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;
       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; // TODO constant?
+      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);
       Fill(n);
@@ -123 +123 @@
       var point = new double[dimension];
       Buffer.BlockCopy(data, sizeof(double) * dimension * newIndex, point, 0, sizeof(double) * dimension);
-      if(!boundary.ContainsPoint(point)) return false;
+      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 +140 @@
       // 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 +165 @@
       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;
@@ -177 +177 @@
 
       // 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]);
-        }
-        index[i] = int.MaxValue;
+        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
       }
 
@@ -192 +192 @@
     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 * 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();
+      for (var i = 0; i < noChildren; i++) correct = correct && children[i].IsCorrect();
       return correct;
     }
@@ -206 +208 @@
     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;
     }
@@ -220 +222 @@
     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 +242 @@
         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 +252 @@
     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];
         }
       }
@@ -269 +270 @@
     #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) {
@@ -275 +276 @@
       dimension = inpData.GetLength(1);
       noChildren = 2;
-      for(uint i = 1; i < dimension; i++) noChildren *= 2;
+      for (uint i = 1; i < dimension; i++) noChildren *= 2;
       data = inpData;
       isLeaf = true;
@@ -316 +317 @@
       }
       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;
       }

branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/TSNEAlgorithm.cs

@@ -285 +285 @@
         } else if((problemData.Dataset as Dataset).VariableHasType<double>(Classes)) {
           var classValues = problemData.Dataset.GetDoubleValues(Classes).ToArray();
-          var max = classValues.Max() + 0.1;     // TODO consts
+          var max = classValues.Max() + 0.1;    
           var min = classValues.Min() - 0.1;
           const int contours = 8;

branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/TSNEStatic.cs

@@ -285 +285 @@
           var minBeta = double.MinValue;
           var maxBeta = double.MaxValue;
-          const double tol = 1e-5;  // TODO: why 1e-5?
+          const double tol = 1e-5;  
 
           // Iterate until we found a good perplexity
           var iter = 0; double sumP = 0;
-          while(!found && iter < 200) {  // TODO 200 iterations always ok?
+          while(!found && iter < 200) {  
 
             // Compute Gaussian kernel row
-            for(var m = 0; m < k; m++) curP[m] = Math.Exp(-beta * distances[m + 1]); // TODO distances m+1?
+            for(var m = 0; m < k; m++) curP[m] = Math.Exp(-beta * distances[m + 1]); 
 
             // Compute entropy of current row
@@ -298 +298 @@
             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]); // TODO: distances m+1?
+            for(var m = 0; m < k; m++) h += beta * (distances[m + 1] * curP[m]); 
             h = h / sumP + Math.Log(sumP);
 

branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/VantagePointTree.cs

@@ -84 +84 @@
     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);
-      Search(root, target, k, heap, double.MaxValue);
+      double 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]);
+        res.Add(items[heap.PeekMinValue().Index]);          // actually max distance
         dist.Add(heap.PeekMinValue().Value);
         heap.RemoveMin();
       }
-      res.Reverse();
+      res.Reverse();  
       dist.Reverse();
       results = res;
@@ -98 +99 @@
     }
 
-    private void Search(Node node, T target, int k, PriorityQueue<double, IndexedItem<double>> heap, double tau = double.MaxValue) {
+    private void Search(Node node, T target, int k, PriorityQueue<double, IndexedItem<double>> heap, ref double tau) {
       if (node == null) return;
       var dist = distance.Get(items[node.index], target);
       if (dist < tau) {
-        if (heap.Size == k) heap.RemoveMin();
-        heap.Insert(-dist, new IndexedItem<double>(node.index, dist));//TODO check if minheap or maxheap should be used here
+        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));      
         if (heap.Size == k) tau = heap.PeekMinValue().Value;
       }
@@ -109 +110 @@
 
       if (dist < node.threshold) {
-        if (dist - tau <= node.threshold) Search(node.left, target, k, heap, tau);   // if there can still be neighbors inside the ball, recursively search left child first
-        if (dist + tau >= node.threshold) Search(node.right, target, k, heap, tau);  // if there can still be neighbors outside the ball, recursively search right child
+        if (dist - tau <= node.threshold) Search(node.left, target, k, heap, ref tau);   // if there can still be neighbors inside the ball, recursively search left child first
+        if (dist + tau >= node.threshold) Search(node.right, target, k, heap, ref tau);  // if there can still be neighbors outside the ball, recursively search right child
       } else {
-        if (dist + tau >= node.threshold) Search(node.right, target, k, heap, tau);  // if there can still be neighbors outside the ball, recursively search right child first
-        if (dist - tau <= node.threshold) Search(node.left, target, k, heap, tau);   // if there can still be neighbors inside the ball, recursively search left child
+        if (dist + tau >= node.threshold) Search(node.right, target, k, heap, ref tau);  // if there can still be neighbors outside the ball, recursively search right child first
+        if (dist - tau <= node.threshold) Search(node.left, target, k, heap, ref tau);   // if there can still be neighbors inside the ball, recursively search left child
       }
     }
@@ -124 +125 @@
       var node = new Node { index = lower };
       var r = new MersenneTwister(); //outer behaviour does not change with the random seed => no need to take the IRandom from the algorithm 
-      if (upper - lower <= 1) return node; // if we did not arrive at leaf yet
+      if (upper - lower <= 1) return node;
 
+      // if we did not arrive at leaf yet
       // Choose an arbitrary point and move it to the start
-      var i = (int)(r.NextDouble() / 1 * (upper - lower - 1)) + lower;
+      var i = (int)(r.NextDouble() * (upper - lower - 1)) + lower;
       items.Swap(lower, i);