[14414] | 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 |
|
---|
| 56 | using System;
|
---|
| 57 | using System.Collections.Generic;
|
---|
| 58 |
|
---|
| 59 | namespace HeuristicLab.Algorithms.DataAnalysis {
|
---|
[14783] | 60 | /// <summary>
|
---|
[14785] | 61 | /// Space partitioning tree (SPTree)
|
---|
[14783] | 62 | /// </summary>
|
---|
[15207] | 63 | internal class SpacePartitioningTree {
|
---|
| 64 | private const uint QtNodeCapacity = 1;
|
---|
[14414] | 65 |
|
---|
[15207] | 66 | #region Fields
|
---|
[14414] | 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
|
---|
[14787] | 73 | private Cell boundary;
|
---|
[14414] | 74 |
|
---|
[14806] | 75 | private double[,] data;
|
---|
[14787] | 76 |
|
---|
[14855] | 77 | // Indices in this space-partitioning tree node, corresponding center-of-mass, and list of all children
|
---|
[14414] | 78 | private double[] centerOfMass;
|
---|
[15207] | 79 | private readonly int[] index = new int[QtNodeCapacity];
|
---|
[14414] | 80 |
|
---|
| 81 | // Children
|
---|
[14785] | 82 | private SpacePartitioningTree[] children;
|
---|
[14414] | 83 | private uint noChildren;
|
---|
[15207] | 84 | #endregion
|
---|
[14414] | 85 |
|
---|
[14785] | 86 | public SpacePartitioningTree(double[,] inpData) {
|
---|
[14414] | 87 | var d = inpData.GetLength(1);
|
---|
| 88 | var n = inpData.GetLength(0);
|
---|
| 89 | var meanY = new double[d];
|
---|
| 90 | var minY = new double[d];
|
---|
[14855] | 91 | for (var i = 0; i < d; i++) minY[i] = double.MaxValue;
|
---|
[14414] | 92 | var maxY = new double[d];
|
---|
[14855] | 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];
|
---|
[14414] | 99 | }
|
---|
| 100 | }
|
---|
[14855] | 101 | for (var i = 0; i < d; i++) meanY[i] /= n;
|
---|
[14414] | 102 | var width = new double[d];
|
---|
[14855] | 103 | for (var i = 0; i < d; i++) width[i] = Math.Max(maxY[i] - meanY[i], meanY[i] - minY[i]) + 1e-5;
|
---|
[15207] | 104 | Init(inpData, meanY, width);
|
---|
[14414] | 105 | Fill(n);
|
---|
| 106 | }
|
---|
| 107 |
|
---|
[15207] | 108 | private SpacePartitioningTree(double[,] inpData, IEnumerable<double> impCorner, IEnumerable<double> impWith) {
|
---|
| 109 | Init(inpData, impCorner, impWith);
|
---|
[14414] | 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];
|
---|
[14806] | 115 | Buffer.BlockCopy(data, sizeof(double) * dimension * newIndex, point, 0, sizeof(double) * dimension);
|
---|
[14855] | 116 | if (!boundary.ContainsPoint(point)) return false;
|
---|
[14414] | 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;
|
---|
[14855] | 121 | for (var i = 0; i < dimension; i++) centerOfMass[i] *= mult1;
|
---|
| 122 | for (var i = 0; i < dimension; i++) centerOfMass[i] += mult2 * point[i];
|
---|
[14414] | 123 |
|
---|
| 124 | // If there is space in this quad tree and it is a leaf, add the object here
|
---|
[15207] | 125 | if (isLeaf && size < QtNodeCapacity) {
|
---|
[14414] | 126 | index[size] = newIndex;
|
---|
| 127 | size++;
|
---|
| 128 | return true;
|
---|
| 129 | }
|
---|
| 130 |
|
---|
[15207] | 131 | // Don't add duplicates
|
---|
[14414] | 132 | var anyDuplicate = false;
|
---|
[14855] | 133 | for (uint n = 0; n < size; n++) {
|
---|
[14414] | 134 | var duplicate = true;
|
---|
[14855] | 135 | for (var d = 0; d < dimension; d++) {
|
---|
| 136 | if (Math.Abs(point[d] - data[index[n], d]) < double.Epsilon) continue;
|
---|
[14414] | 137 | duplicate = false; break;
|
---|
| 138 | }
|
---|
| 139 | anyDuplicate = anyDuplicate | duplicate;
|
---|
| 140 | }
|
---|
[14855] | 141 | if (anyDuplicate) return true;
|
---|
[14414] | 142 |
|
---|
| 143 | // Otherwise, we need to subdivide the current cell
|
---|
[14855] | 144 | if (isLeaf) Subdivide();
|
---|
[14414] | 145 | // Find out where the point can be inserted
|
---|
[14855] | 146 | for (var i = 0; i < noChildren; i++) {
|
---|
| 147 | if (children[i].Insert(newIndex)) return true;
|
---|
[14414] | 148 | }
|
---|
| 149 |
|
---|
| 150 | // Otherwise, the point cannot be inserted (this should never happen)
|
---|
| 151 | return false;
|
---|
| 152 | }
|
---|
| 153 |
|
---|
[14806] | 154 | public void ComputeNonEdgeForces(int pointIndex, double theta, double[] negF, ref double sumQ) {
|
---|
[14414] | 155 | // Make sure that we spend no time on empty nodes or self-interactions
|
---|
[14855] | 156 | if (cumulativeSize == 0 || (isLeaf && size == 1 && index[0] == pointIndex)) return;
|
---|
[14414] | 157 |
|
---|
| 158 | // Compute distance between point and center-of-mass
|
---|
| 159 | var D = .0;
|
---|
[15207] | 160 | var buff = new double[dimension];
|
---|
[14855] | 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];
|
---|
[14414] | 163 |
|
---|
| 164 | // Check whether we can use this node as a "summary"
|
---|
| 165 | var maxWidth = 0.0;
|
---|
[14855] | 166 | for (var d = 0; d < dimension; d++) {
|
---|
[14414] | 167 | var curWidth = boundary.GetWidth(d);
|
---|
[15207] | 168 | maxWidth = maxWidth > curWidth ? maxWidth : curWidth;
|
---|
[14414] | 169 | }
|
---|
[14855] | 170 | if (isLeaf || maxWidth / Math.Sqrt(D) < theta) {
|
---|
[14414] | 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;
|
---|
[14788] | 175 | sumQ += mult;
|
---|
[14414] | 176 | mult *= D;
|
---|
[14855] | 177 | for (var d = 0; d < dimension; d++) negF[d] += mult * buff[d];
|
---|
[14414] | 178 | } else {
|
---|
| 179 |
|
---|
| 180 | // Recursively apply Barnes-Hut to children
|
---|
[14855] | 181 | for (var i = 0; i < noChildren; i++) children[i].ComputeNonEdgeForces(pointIndex, theta, negF, ref sumQ);
|
---|
[14414] | 182 | }
|
---|
| 183 | }
|
---|
| 184 |
|
---|
[15207] | 185 | public static void ComputeEdgeForces(int[] rowP, int[] colP, double[] valP, int n, double[,] posF, double[,] data, int dimension) {
|
---|
[14414] | 186 | // Loop over all edges in the graph
|
---|
[14855] | 187 | for (var k = 0; k < n; k++) {
|
---|
| 188 | for (var i = rowP[k]; i < rowP[k + 1]; i++) {
|
---|
[14414] | 189 |
|
---|
| 190 | // Compute pairwise distance and Q-value
|
---|
[14806] | 191 | // uses squared distance
|
---|
[14414] | 192 | var d = 1.0;
|
---|
[15207] | 193 | var buff = new double[dimension];
|
---|
[14855] | 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];
|
---|
[14414] | 196 | d = valP[i] / d;
|
---|
| 197 |
|
---|
| 198 | // Sum positive force
|
---|
[14855] | 199 | for (var j = 0; j < dimension; j++) posF[k, j] += d * buff[j];
|
---|
[14414] | 200 | }
|
---|
| 201 | }
|
---|
| 202 | }
|
---|
| 203 |
|
---|
| 204 | #region Helpers
|
---|
| 205 | private void Fill(int n) {
|
---|
[14855] | 206 | for (var i = 0; i < n; i++) Insert(i);
|
---|
[14414] | 207 | }
|
---|
[15207] | 208 |
|
---|
| 209 | private void Init(double[,] inpData, IEnumerable<double> inpCorner, IEnumerable<double> inpWidth) {
|
---|
[14414] | 210 | dimension = inpData.GetLength(1);
|
---|
| 211 | noChildren = 2;
|
---|
[14855] | 212 | for (uint i = 1; i < dimension; i++) noChildren *= 2;
|
---|
[14806] | 213 | data = inpData;
|
---|
[14414] | 214 | isLeaf = true;
|
---|
| 215 | size = 0;
|
---|
| 216 | cumulativeSize = 0;
|
---|
| 217 | boundary = new Cell((uint)dimension);
|
---|
[15207] | 218 |
|
---|
[14414] | 219 | inpCorner.ForEach((i, x) => boundary.SetCorner(i, x));
|
---|
| 220 | inpWidth.ForEach((i, x) => boundary.SetWidth(i, x));
|
---|
| 221 |
|
---|
[14785] | 222 | children = new SpacePartitioningTree[noChildren];
|
---|
[14414] | 223 | centerOfMass = new double[dimension];
|
---|
[15207] | 224 | }
|
---|
[14414] | 225 |
|
---|
[15207] | 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;
|
---|
[14414] | 252 | }
|
---|
| 253 | #endregion
|
---|
| 254 |
|
---|
[14806] | 255 | private class Cell {
|
---|
[14787] | 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) {
|
---|
[14855] | 279 | for (var d = 0; d < dimension; d++)
|
---|
| 280 | if (corner[d] - width[d] > point[d] || corner[d] + width[d] < point[d]) return false;
|
---|
[14787] | 281 | return true;
|
---|
| 282 | }
|
---|
| 283 | }
|
---|
[14414] | 284 | }
|
---|
| 285 | }
|
---|