Free cookie consent management tool by TermsFeed Policy Generator

source: branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/SPtree.cs @ 14518

Last change on this file since 14518 was 14518, checked in by bwerth, 7 years ago

#2700 TSNEAnalysis is now a BasicAlg, hid some Parameters, added optional data normalization to make TSNE scaling-invariant

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