source: branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/SpacePartitioningTree.cs @ 14855

Last change on this file since 14855 was 14855, checked in by gkronber, 2 years ago

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

File size: 12.2 KB
Line 
1#region License Information
2/* HeuristicLab
3 * Copyright (C) 2002-2016 Heuristic and Evolutionary Algorithms Laboratory (HEAL)
4 *
5 * This file is part of HeuristicLab.
6 *
7 * HeuristicLab is free software: you can redistribute it and/or modify
8 * it under the terms of the GNU General Public License as published by
9 * the Free Software Foundation, either version 3 of the License, or
10 * (at your option) any later version.
11 *
12 * HeuristicLab is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
15 * GNU General Public License for more details.
16 *
17 * You should have received a copy of the GNU General Public License
18 * along with HeuristicLab. If not, see <http://www.gnu.org/licenses/>.
19 */
20
21//Code is based on an implementation from Laurens van der Maaten
22
23/*
24*
25* Copyright (c) 2014, Laurens van der Maaten (Delft University of Technology)
26* All rights reserved.
27*
28* Redistribution and use in source and binary forms, with or without
29* modification, are permitted provided that the following conditions are met:
30* 1. Redistributions of source code must retain the above copyright
31*    notice, this list of conditions and the following disclaimer.
32* 2. Redistributions in binary form must reproduce the above copyright
33*    notice, this list of conditions and the following disclaimer in the
34*    documentation and/or other materials provided with the distribution.
35* 3. All advertising materials mentioning features or use of this software
36*    must display the following acknowledgement:
37*    This product includes software developed by the Delft University of Technology.
38* 4. Neither the name of the Delft University of Technology nor the names of
39*    its contributors may be used to endorse or promote products derived from
40*    this software without specific prior written permission.
41*
42* THIS SOFTWARE IS PROVIDED BY LAURENS VAN DER MAATEN ''AS IS'' AND ANY EXPRESS
43* OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
44* OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO
45* EVENT SHALL LAURENS VAN DER MAATEN BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
46* SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
47* PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR
48* BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
49* CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING
50* IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY
51* OF SUCH DAMAGE.
52*
53*/
54#endregion
55
56using System;
57using System.Collections.Generic;
58using System.Linq;
59using HeuristicLab.Common;
60
61namespace HeuristicLab.Algorithms.DataAnalysis {
62  /// <summary>
63  /// Space partitioning tree (SPTree)
64  /// </summary>
65  public class SpacePartitioningTree : ISpacePartitioningTree {
66    private const uint QT_NODE_CAPACITY = 1;
67
68    private double[] buff;
69    private SpacePartitioningTree parent;
70    private int dimension;
71    private bool isLeaf;
72    private uint size;
73    private uint cumulativeSize;
74
75    // Axis-aligned bounding box stored as a center with half-dimensions to represent the boundaries of this quad tree
76    private Cell boundary;
77
78    private double[,] data;
79
80    // Indices in this space-partitioning tree node, corresponding center-of-mass, and list of all children
81    private double[] centerOfMass;
82    private readonly int[] index = new int[QT_NODE_CAPACITY];
83
84    // Children
85    private SpacePartitioningTree[] children;
86    private uint noChildren;
87
88    public SpacePartitioningTree(double[,] inpData) {
89      var d = inpData.GetLength(1);
90      var n = inpData.GetLength(0);
91      var meanY = new double[d];
92      var minY = new double[d];
93      for (var i = 0; i < d; i++) minY[i] = double.MaxValue;
94      var maxY = new double[d];
95      for (var i = 0; i < d; i++) maxY[i] = double.MinValue;
96      for (uint i = 0; i < n; i++) {
97        for (uint j = 0; j < d; j++) {
98          meanY[j] += inpData[i, j];
99          if (inpData[i, j] < minY[j]) minY[j] = inpData[i, j];
100          if (inpData[i, j] > maxY[j]) maxY[j] = inpData[i, j];
101        }
102      }
103      for (var i = 0; i < d; i++) meanY[i] /= n;
104      var width = new double[d];
105      for (var i = 0; i < d; i++) width[i] = Math.Max(maxY[i] - meanY[i], meanY[i] - minY[i]) + 1e-5;
106      Init(null, inpData, meanY, width);
107      Fill(n);
108    }
109
110    public SpacePartitioningTree(double[,] inpData, IEnumerable<double> impCorner, IEnumerable<double> impWith) {
111      Init(null, inpData, impCorner, impWith);
112    }
113    public SpacePartitioningTree(SpacePartitioningTree parent, double[,] inpData, IEnumerable<double> impCorner, IEnumerable<double> impWith) {
114      Init(parent, inpData, impCorner, impWith);
115    }
116
117    public ISpacePartitioningTree GetParent() {
118      return parent;
119    }
120
121    public bool Insert(int newIndex) {
122      // Ignore objects which do not belong in this quad tree
123      var point = new double[dimension];
124      Buffer.BlockCopy(data, sizeof(double) * dimension * newIndex, point, 0, sizeof(double) * dimension);
125      if (!boundary.ContainsPoint(point)) return false;
126      cumulativeSize++;
127      // Online update of cumulative size and center-of-mass
128      var mult1 = (double)(cumulativeSize - 1) / cumulativeSize;
129      var mult2 = 1.0 / cumulativeSize;
130      for (var i = 0; i < dimension; i++) centerOfMass[i] *= mult1;
131      for (var i = 0; i < dimension; i++) centerOfMass[i] += mult2 * point[i];
132
133      // If there is space in this quad tree and it is a leaf, add the object here
134      if (isLeaf && size < QT_NODE_CAPACITY) {
135        index[size] = newIndex;
136        size++;
137        return true;
138      }
139
140      // Don't add duplicates for now (this is not very nice)
141      var anyDuplicate = false;
142      for (uint n = 0; n < size; n++) {
143        var duplicate = true;
144        for (var d = 0; d < dimension; d++) {
145          if (Math.Abs(point[d] - data[index[n], d]) < double.Epsilon) continue;
146          duplicate = false; break;
147        }
148        anyDuplicate = anyDuplicate | duplicate;
149      }
150      if (anyDuplicate) return true;
151
152      // Otherwise, we need to subdivide the current cell
153      if (isLeaf) Subdivide();
154      // Find out where the point can be inserted
155      for (var i = 0; i < noChildren; i++) {
156        if (children[i].Insert(newIndex)) return true;
157      }
158
159      // Otherwise, the point cannot be inserted (this should never happen)
160      return false;
161    }
162
163    public void Subdivide() {
164      // Create new children
165      var newCorner = new double[dimension];
166      var newWidth = new double[dimension];
167      for (var i = 0; i < noChildren; i++) {
168        var div = 1;
169        for (var d = 0; d < dimension; d++) {
170          newWidth[d] = .5 * boundary.GetWidth(d);
171          if ((i / div) % 2 == 1) newCorner[d] = boundary.GetCorner(d) - .5 * boundary.GetWidth(d);
172          else newCorner[d] = boundary.GetCorner(d) + .5 * boundary.GetWidth(d);
173          div *= 2;
174        }
175        children[i] = new SpacePartitioningTree(this, data, newCorner, newWidth);
176      }
177
178      // Move existing points to correct children
179      for (var i = 0; i < size; i++) {
180        var success = false;
181        for (var j = 0; j < noChildren; j++) {
182          if (!success) success = children[j].Insert(index[i]);
183        }
184        index[i] = -1; // as in tSNE implementation by van der Maaten
185      }
186
187      // Empty parent node
188      size = 0;
189      isLeaf = false;
190    }
191
192    public bool IsCorrect() {
193      var row = new double[dimension];
194      for (var n = 0; n < size; n++) {
195        Buffer.BlockCopy(data, sizeof(double) * dimension * index[n], row, 0, sizeof(double) * dimension);
196        if (!boundary.ContainsPoint(row)) return false;
197      }
198      if (isLeaf) return true;
199      var correct = true;
200      for (var i = 0; i < noChildren; i++) correct = correct && children[i].IsCorrect();
201      return correct;
202    }
203
204    public void GetAllIndices(int[] indices) {
205      GetAllIndices(indices, 0);
206    }
207
208    public int GetAllIndices(int[] indices, int loc) {
209      // Gather indices in current quadrant
210      for (var i = 0; i < size; i++) indices[loc + i] = index[i];
211      loc += (int)size;
212      // Gather indices in children
213      if (isLeaf) return loc;
214      for (var i = 0; i < noChildren; i++) loc = children[i].GetAllIndices(indices, loc);
215      return loc;
216    }
217
218    public int GetDepth() {
219      return isLeaf ? 1 : 1 + children.Max(x => x.GetDepth());
220    }
221
222    public void ComputeNonEdgeForces(int pointIndex, double theta, double[] negF, ref double sumQ) {
223      // Make sure that we spend no time on empty nodes or self-interactions
224      if (cumulativeSize == 0 || (isLeaf && size == 1 && index[0] == pointIndex)) return;
225
226      // Compute distance between point and center-of-mass
227      var D = .0;
228      for (var d = 0; d < dimension; d++) buff[d] = data[pointIndex, d] - centerOfMass[d];
229      for (var d = 0; d < dimension; d++) D += buff[d] * buff[d];
230
231      // Check whether we can use this node as a "summary"
232      var maxWidth = 0.0;
233      for (var d = 0; d < dimension; d++) {
234        var curWidth = boundary.GetWidth(d);
235        maxWidth = (maxWidth > curWidth) ? maxWidth : curWidth;
236      }
237      if (isLeaf || maxWidth / Math.Sqrt(D) < theta) {
238
239        // Compute and add t-SNE force between point and current node
240        D = 1.0 / (1.0 + D);
241        var mult = cumulativeSize * D;
242        sumQ += mult;
243        mult *= D;
244        for (var d = 0; d < dimension; d++) negF[d] += mult * buff[d];
245      } else {
246
247        // Recursively apply Barnes-Hut to children
248        for (var i = 0; i < noChildren; i++) children[i].ComputeNonEdgeForces(pointIndex, theta, negF, ref sumQ);
249      }
250    }
251
252    public void ComputeEdgeForces(int[] rowP, int[] colP, double[] valP, int n, double[,] posF) {
253      // Loop over all edges in the graph
254      for (var k = 0; k < n; k++) {
255        for (var i = rowP[k]; i < rowP[k + 1]; i++) {
256
257          // Compute pairwise distance and Q-value
258          // uses squared distance
259          var d = 1.0;
260          for (var j = 0; j < dimension; j++) buff[j] = data[k, j] - data[colP[i], j];
261          for (var j = 0; j < dimension; j++) d += buff[j] * buff[j];
262          d = valP[i] / d;
263
264          // Sum positive force
265          for (var j = 0; j < dimension; j++) posF[k, j] += d * buff[j];
266        }
267      }
268    }
269
270    #region Helpers
271    private void Fill(int n) {
272      for (var i = 0; i < n; i++) Insert(i);
273    }
274    private void Init(SpacePartitioningTree p, double[,] inpData, IEnumerable<double> inpCorner, IEnumerable<double> inpWidth) {
275      parent = p;
276      dimension = inpData.GetLength(1);
277      noChildren = 2;
278      for (uint i = 1; i < dimension; i++) noChildren *= 2;
279      data = inpData;
280      isLeaf = true;
281      size = 0;
282      cumulativeSize = 0;
283      boundary = new Cell((uint)dimension);
284      inpCorner.ForEach((i, x) => boundary.SetCorner(i, x));
285      inpWidth.ForEach((i, x) => boundary.SetWidth(i, x));
286
287      children = new SpacePartitioningTree[noChildren];
288      centerOfMass = new double[dimension];
289      buff = new double[dimension];
290
291    }
292    #endregion
293
294
295    private class Cell {
296      private readonly uint dimension;
297      private readonly double[] corner;
298      private readonly double[] width;
299
300      public Cell(uint inpDimension) {
301        dimension = inpDimension;
302        corner = new double[dimension];
303        width = new double[dimension];
304      }
305
306      public double GetCorner(int d) {
307        return corner[d];
308      }
309      public double GetWidth(int d) {
310        return width[d];
311      }
312      public void SetCorner(int d, double val) {
313        corner[d] = val;
314      }
315      public void SetWidth(int d, double val) {
316        width[d] = val;
317      }
318      public bool ContainsPoint(double[] point) {
319        for (var d = 0; d < dimension; d++)
320          if (corner[d] - width[d] > point[d] || corner[d] + width[d] < point[d]) return false;
321        return true;
322      }
323    }
324  }
325}
Note: See TracBrowser for help on using the repository browser.