1  #region License Information


2  /* HeuristicLab


3  * Copyright (C) 20022016 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  using System.Linq;


59  using HeuristicLab.Common;


60 


61  namespace 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  // Axisaligned bounding box stored as a center with halfdimensions to represent the boundaries of this quad tree


76  private Cell boundary;


77 


78  // Indices in this spacepartitioning tree node, corresponding centerofmass, and list of all children


79  private double[,] data;


80 


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]  maxY[i]) + 1e5; // TODO constant?


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 centerofmass


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] = int.MaxValue;


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++) Buffer.BlockCopy(data, sizeof(double) * dimension * n, row, 0, sizeof(double) * dimension);


195  if(!boundary.ContainsPoint(row)) return false;


196  if(isLeaf) return true;


197  var correct = true;


198  for(var i = 0; i < noChildren; i++) correct = correct && children[i].IsCorrect();


199  return correct;


200  }


201 


202  public void GetAllIndices(int[] indices) {


203  GetAllIndices(indices, 0);


204  }


205 


206  public int GetAllIndices(int[] indices, int loc) {


207  // Gather indices in current quadrant


208  for(var i = 0; i < size; i++) indices[loc + i] = index[i];


209  loc += (int)size;


210  // Gather indices in children


211  if(isLeaf) return loc;


212  for(var i = 0; i < noChildren; i++) loc = children[i].GetAllIndices(indices, loc);


213  return loc;


214  }


215 


216  public int GetDepth() {


217  return isLeaf ? 1 : 1 + children.Max(x => x.GetDepth());


218  }


219 


220  public void ComputeNonEdgeForces(int pointIndex, double theta, double[] negF, ref double sumQ) {


221  // Make sure that we spend no time on empty nodes or selfinteractions


222  if(cumulativeSize == 0  (isLeaf && size == 1 && index[0] == pointIndex)) return;


223 


224  // Compute distance between point and centerofmass


225  // TODO: squared distance with normalized axes is used here!


226  var D = .0;


227  for(var d = 0; d < dimension; d++) buff[d] = data[pointIndex, d]  centerOfMass[d];


228  for(var d = 0; d < dimension; d++) D += buff[d] * buff[d];


229 


230  // Check whether we can use this node as a "summary"


231  var maxWidth = 0.0;


232  for(var d = 0; d < dimension; d++) {


233  var curWidth = boundary.GetWidth(d);


234  maxWidth = (maxWidth > curWidth) ? maxWidth : curWidth;


235  }


236  if(isLeaf  maxWidth / Math.Sqrt(D) < theta) {


237 


238  // Compute and add tSNE force between point and current node


239  D = 1.0 / (1.0 + D);


240  var mult = cumulativeSize * D;


241  sumQ += mult;


242  mult *= D;


243  for(var d = 0; d < dimension; d++) negF[d] += mult * buff[d];


244  } else {


245 


246  // Recursively apply BarnesHut to children


247  for(var i = 0; i < noChildren; i++) children[i].ComputeNonEdgeForces(pointIndex, theta, negF, ref sumQ);


248  }


249  }


250 


251  public void ComputeEdgeForces(int[] rowP, int[] colP, double[] valP, int n, double[,] posF) {


252  // Loop over all edges in the graph


253  for(var k = 0; k < n; k++) {


254  for(var i = rowP[k]; i < rowP[k + 1]; i++) {


255 


256  // Compute pairwise distance and Qvalue


257  // uses squared distance


258  var d = 1.0;


259  for(var j = 0; j < dimension; j++) buff[j] = data[k, j]  data[colP[i], j];


260  for(var j = 0; j < dimension; j++) d += buff[j] * buff[j];


261  d = valP[i] / d;


262 


263  // Sum positive force


264  for(var j = 0; j < dimension; j++) posF[k, j] += d * buff[j];


265  }


266  }


267  }


268 


269  #region Helpers


270  private void Fill(int n) {


271  for(var i = 0; i < n; i++) Insert(i);


272  }


273  private void Init(SpacePartitioningTree p, double[,] inpData, IEnumerable<double> inpCorner, IEnumerable<double> inpWidth) {


274  parent = p;


275  dimension = inpData.GetLength(1);


276  noChildren = 2;


277  for(uint i = 1; i < dimension; i++) noChildren *= 2;


278  data = inpData;


279  isLeaf = true;


280  size = 0;


281  cumulativeSize = 0;


282  boundary = new Cell((uint)dimension);


283  inpCorner.ForEach((i, x) => boundary.SetCorner(i, x));


284  inpWidth.ForEach((i, x) => boundary.SetWidth(i, x));


285 


286  children = new SpacePartitioningTree[noChildren];


287  centerOfMass = new double[dimension];


288  buff = new double[dimension];


289 


290  }


291  #endregion


292 


293 


294  private class Cell {


295  private readonly uint dimension;


296  private readonly double[] corner;


297  private readonly double[] width;


298 


299  public Cell(uint inpDimension) {


300  dimension = inpDimension;


301  corner = new double[dimension];


302  width = new double[dimension];


303  }


304 


305  public double GetCorner(int d) {


306  return corner[d];


307  }


308  public double GetWidth(int d) {


309  return width[d];


310  }


311  public void SetCorner(int d, double val) {


312  corner[d] = val;


313  }


314  public void SetWidth(int d, double val) {


315  width[d] = val;


316  }


317  public bool ContainsPoint(double[] point) {


318  for(var d = 0; d < dimension; d++)


319  if(corner[d]  width[d] > point[d]  corner[d] + width[d] < point[d]) return false;


320  return true;


321  }


322  }


323  }


324  }

