source: branches/2994-AutoDiffForIntervals/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/SpacePartitioningTree.cs @ 17209

Last change on this file since 17209 was 17209, checked in by gkronber, 6 weeks ago

#2994: merged r17132:17198 from trunk to branch

File size: 10.8 KB
Line 
1#region License Information
2/* HeuristicLab
3 * Copyright (C) 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;
58
59namespace HeuristicLab.Algorithms.DataAnalysis {
60  /// <summary>
61  /// Space partitioning tree (SPTree)
62  /// </summary>
63  internal class SpacePartitioningTree {
64    private const uint QtNodeCapacity = 1;
65
66    #region Fields
67    private int dimension;
68    private bool isLeaf;
69    private uint size;
70    private uint cumulativeSize;
71
72    // Axis-aligned bounding box stored as a center with half-dimensions to represent the boundaries of this quad tree
73    private Cell boundary;
74
75    private double[,] data;
76
77    // Indices in this space-partitioning tree node, corresponding center-of-mass, and list of all children
78    private double[] centerOfMass;
79    private readonly int[] index = new int[QtNodeCapacity];
80
81    // Children
82    private SpacePartitioningTree[] children;
83    private uint noChildren;
84    #endregion
85
86    public SpacePartitioningTree(double[,] inpData) {
87      var d = inpData.GetLength(1);
88      var n = inpData.GetLength(0);
89      var meanY = new double[d];
90      var minY = new double[d];
91      for (var i = 0; i < d; i++) minY[i] = double.MaxValue;
92      var maxY = new double[d];
93      for (var i = 0; i < d; i++) maxY[i] = double.MinValue;
94      for (uint i = 0; i < n; i++) {
95        for (uint j = 0; j < d; j++) {
96          meanY[j] += inpData[i, j];
97          if (inpData[i, j] < minY[j]) minY[j] = inpData[i, j];
98          if (inpData[i, j] > maxY[j]) maxY[j] = inpData[i, j];
99        }
100      }
101      for (var i = 0; i < d; i++) meanY[i] /= n;
102      var width = new double[d];
103      for (var i = 0; i < d; i++) width[i] = Math.Max(maxY[i] - meanY[i], meanY[i] - minY[i]) + 1e-5;
104      Init(inpData, meanY, width);
105      Fill(n);
106    }
107
108    private SpacePartitioningTree(double[,] inpData, IEnumerable<double> impCorner, IEnumerable<double> impWith) {
109      Init(inpData, impCorner, impWith);
110    }
111
112    public bool Insert(int newIndex) {
113      // Ignore objects which do not belong in this quad tree
114      var point = new double[dimension];
115      Buffer.BlockCopy(data, sizeof(double) * dimension * newIndex, point, 0, sizeof(double) * dimension);
116      if (!boundary.ContainsPoint(point)) return false;
117      cumulativeSize++;
118      // Online update of cumulative size and center-of-mass
119      var mult1 = (double)(cumulativeSize - 1) / cumulativeSize;
120      var mult2 = 1.0 / cumulativeSize;
121      for (var i = 0; i < dimension; i++) centerOfMass[i] *= mult1;
122      for (var i = 0; i < dimension; i++) centerOfMass[i] += mult2 * point[i];
123
124      // If there is space in this quad tree and it is a leaf, add the object here
125      if (isLeaf && size < QtNodeCapacity) {
126        index[size] = newIndex;
127        size++;
128        return true;
129      }
130
131      // Don't add duplicates
132      var anyDuplicate = false;
133      for (uint n = 0; n < size; n++) {
134        var duplicate = true;
135        for (var d = 0; d < dimension; d++) {
136          if (Math.Abs(point[d] - data[index[n], d]) < double.Epsilon) continue;
137          duplicate = false; break;
138        }
139        anyDuplicate = anyDuplicate | duplicate;
140      }
141      if (anyDuplicate) return true;
142
143      // Otherwise, we need to subdivide the current cell
144      if (isLeaf) Subdivide();
145      // Find out where the point can be inserted
146      for (var i = 0; i < noChildren; i++) {
147        if (children[i].Insert(newIndex)) return true;
148      }
149
150      // Otherwise, the point cannot be inserted (this should never happen)
151      return false;
152    }
153
154    public void ComputeNonEdgeForces(int pointIndex, double theta, double[] negF, ref double sumQ) {
155      // Make sure that we spend no time on empty nodes or self-interactions
156      if (cumulativeSize == 0 || (isLeaf && size == 1 && index[0] == pointIndex)) return;
157
158      // Compute distance between point and center-of-mass
159      var D = .0;
160      var buff = new double[dimension];
161      for (var d = 0; d < dimension; d++) buff[d] = data[pointIndex, d] - centerOfMass[d];
162      for (var d = 0; d < dimension; d++) D += buff[d] * buff[d];
163
164      // Check whether we can use this node as a "summary"
165      var maxWidth = 0.0;
166      for (var d = 0; d < dimension; d++) {
167        var curWidth = boundary.GetWidth(d);
168        maxWidth = maxWidth > curWidth ? maxWidth : curWidth;
169      }
170      if (isLeaf || maxWidth / Math.Sqrt(D) < theta) {
171
172        // Compute and add t-SNE force between point and current node
173        D = 1.0 / (1.0 + D);
174        var mult = cumulativeSize * D;
175        sumQ += mult;
176        mult *= D;
177        for (var d = 0; d < dimension; d++) negF[d] += mult * buff[d];
178      } else {
179
180        // Recursively apply Barnes-Hut to children
181        for (var i = 0; i < noChildren; i++) children[i].ComputeNonEdgeForces(pointIndex, theta, negF, ref sumQ);
182      }
183    }
184
185    public static void ComputeEdgeForces(int[] rowP, int[] colP, double[] valP, int n, double[,] posF, double[,] data, int dimension) {
186      // Loop over all edges in the graph
187      for (var k = 0; k < n; k++) {
188        for (var i = rowP[k]; i < rowP[k + 1]; i++) {
189
190          // Compute pairwise distance and Q-value
191          // uses squared distance
192          var d = 1.0;
193          var buff = new double[dimension];
194          for (var j = 0; j < dimension; j++) buff[j] = data[k, j] - data[colP[i], j];
195          for (var j = 0; j < dimension; j++) d += buff[j] * buff[j];
196          d = valP[i] / d;
197
198          // Sum positive force
199          for (var j = 0; j < dimension; j++) posF[k, j] += d * buff[j];
200        }
201      }
202    }
203
204    #region Helpers
205    private void Fill(int n) {
206      for (var i = 0; i < n; i++) Insert(i);
207    }
208
209    private void Init(double[,] inpData, IEnumerable<double> inpCorner, IEnumerable<double> inpWidth) {
210      dimension = inpData.GetLength(1);
211      noChildren = 2;
212      for (uint i = 1; i < dimension; i++) noChildren *= 2;
213      data = inpData;
214      isLeaf = true;
215      size = 0;
216      cumulativeSize = 0;
217      boundary = new Cell((uint)dimension);
218
219      inpCorner.ForEach((i, x) => boundary.SetCorner(i, x));
220      inpWidth.ForEach((i, x) => boundary.SetWidth(i, x));
221
222      children = new SpacePartitioningTree[noChildren];
223      centerOfMass = new double[dimension];
224    }
225
226    private void Subdivide() {
227      // Create new children
228      var newCorner = new double[dimension];
229      var newWidth = new double[dimension];
230      for (var i = 0; i < noChildren; i++) {
231        var div = 1;
232        for (var d = 0; d < dimension; d++) {
233          newWidth[d] = .5 * boundary.GetWidth(d);
234          if (i / div % 2 == 1) newCorner[d] = boundary.GetCorner(d) - .5 * boundary.GetWidth(d);
235          else newCorner[d] = boundary.GetCorner(d) + .5 * boundary.GetWidth(d);
236          div *= 2;
237        }
238        children[i] = new SpacePartitioningTree(data, newCorner, newWidth);
239      }
240
241      // Move existing points to correct children
242      for (var i = 0; i < size; i++) {
243        var success = false;
244        for (var j = 0; j < noChildren; j++) {
245          if (!success) success = children[j].Insert(index[i]);
246        }
247        index[i] = -1; // as in tSNE implementation by van der Maaten
248      }
249      // Empty parent node
250      size = 0;
251      isLeaf = false;
252    }
253    #endregion
254
255    private class Cell {
256      private readonly uint dimension;
257      private readonly double[] corner;
258      private readonly double[] width;
259
260      public Cell(uint inpDimension) {
261        dimension = inpDimension;
262        corner = new double[dimension];
263        width = new double[dimension];
264      }
265
266      public double GetCorner(int d) {
267        return corner[d];
268      }
269      public double GetWidth(int d) {
270        return width[d];
271      }
272      public void SetCorner(int d, double val) {
273        corner[d] = val;
274      }
275      public void SetWidth(int d, double val) {
276        width[d] = val;
277      }
278      public bool ContainsPoint(double[] point) {
279        for (var d = 0; d < dimension; d++)
280          if (corner[d] - width[d] > point[d] || corner[d] + width[d] < point[d]) return false;
281        return true;
282      }
283    }
284  }
285}
Note: See TracBrowser for help on using the repository browser.