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

#region License Information
/* HeuristicLab
 * Copyright (C) 2002-2016 Heuristic and Evolutionary Algorithms Laboratory (HEAL)
 *
 * This file is part of HeuristicLab.
 *
 * HeuristicLab is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * HeuristicLab is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with HeuristicLab. If not, see <http://www.gnu.org/licenses/>.
 */

//Code is based on an implementation from Laurens van der Maaten

/*
*
* Copyright (c) 2014, Laurens van der Maaten (Delft University of Technology)
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions are met:
* 1. Redistributions of source code must retain the above copyright
*    notice, this list of conditions and the following disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright
*    notice, this list of conditions and the following disclaimer in the
*    documentation and/or other materials provided with the distribution.
* 3. All advertising materials mentioning features or use of this software
*    must display the following acknowledgement:
*    This product includes software developed by the Delft University of Technology.
* 4. Neither the name of the Delft University of Technology nor the names of
*    its contributors may be used to endorse or promote products derived from
*    this software without specific prior written permission.
*
* THIS SOFTWARE IS PROVIDED BY LAURENS VAN DER MAATEN ''AS IS'' AND ANY EXPRESS
* OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
* OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO
* EVENT SHALL LAURENS VAN DER MAATEN BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
* SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
* PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR
* BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
* CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING
* IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY
* OF SUCH DAMAGE.
*
*/
#endregion

using System;
using System.Collections.Generic;
using System.Linq;
using HeuristicLab.Common;
using HeuristicLab.Core;
using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
using HeuristicLab.Random;

namespace HeuristicLab.Algorithms.DataAnalysis {
  [StorableClass]
  public class VPTree<T> : DeepCloneable, IVPTree<T> where T : class, IDeepCloneable {
    #region properties
    [Storable]
    private List<T> items;
    [Storable]
    private double tau;
    [Storable]
    private Node root;
    [Storable]
    private IDistance<T> distance;
    #endregion

    #region HLConstructors & Cloning
    [StorableConstructor]
    protected VPTree(bool deserializing) { }
    protected VPTree(VPTree<T> original, Cloner cloner)
      : base(original, cloner) {
      items = original.items.Select(cloner.Clone).ToList();
      tau = original.tau;
      root = cloner.Clone(original.root);
      distance = cloner.Clone(distance);
    }
    public override IDeepCloneable Clone(Cloner cloner) { return new VPTree<T>(this, cloner); }
    public VPTree(IDistance<T> distance) {
      root = null;
      this.distance = distance;
    }
    #endregion

    public void Create(IEnumerable<T> items) {
      this.items = items.Select(x => x).ToList();
      root = BuildFromPoints(0, this.items.Count);
    }

    public void Search(T target, int k, out List<T> results, out List<double> distances) {
      IHeap<double, HeapItem> heap = new PriorityQueue<double, HeapItem>(double.MaxValue, double.MinValue, k);
      tau = double.MaxValue;
      Search(root, target, k, heap);
      results = new List<T>();
      distances = new List<double>();
      while (heap.Size > 0) {
        results.Add(items[heap.PeekMinValue().Index]);
        distances.Add(heap.PeekMinValue().Dist);
        heap.RemoveMin();
      }
      results.Reverse();
      distances.Reverse();
    }

    private Node BuildFromPoints(int lower, int upper) {
      if (upper == lower)      // indicates that we're done here!
        return null;

      // Lower index is center of current node
      var node = new Node { index = lower };
      var r = new MersenneTwister(); //outer behaviour does not change with the random seed => no need to take the IRandom from the algorithm 
      if (upper - lower <= 1) return node; // if we did not arrive at leaf yet

      // Choose an arbitrary point and move it to the start
      var i = (int)(r.NextDouble() / 1 * (upper - lower - 1)) + lower;
      items.Swap(lower, i);

      // Partition around the median distance
      var median = (upper + lower) / 2;
      items.NthElement(lower + 1, upper - 1, median, distance.GetDistanceComparer(items[lower]));

      // Threshold of the new node will be the distance to the median
      node.threshold = distance.Get(items[lower], items[median]);

      // Recursively build tree
      node.index = lower;
      node.left = BuildFromPoints(lower + 1, median);
      node.right = BuildFromPoints(median, upper);

      // Return result
      return node;
    }

    private void Search(Node node, T target, int k, IHeap<double, HeapItem> heap) {
      if (node == null) return;
      var dist = distance.Get(items[node.index], target);
      if (dist < tau) {
        if (heap.Size == k) heap.RemoveMin();
        heap.Insert(-dist, new HeapItem(node.index, dist));//TODO check if minheap or maxheap schould be used here
        if (heap.Size == k) tau = heap.PeekMinValue().Dist;
      }
      if (node.left == null && node.right == null) return;

      if (dist < node.threshold) {
        if (dist - tau <= node.threshold) Search(node.left, target, k, heap);   // if there can still be neighbors inside the ball, recursively search left child first
        if (dist + tau >= node.threshold) Search(node.right, target, k, heap);  // if there can still be neighbors outside the ball, recursively search right child
      } else {
        if (dist + tau >= node.threshold) Search(node.right, target, k, heap);  // if there can still be neighbors outside the ball, recursively search right child first
        if (dist - tau <= node.threshold) Search(node.left, target, k, heap);   // if there can still be neighbors inside the ball, recursively search left child
      }

    }

    [StorableClass]
    private sealed class Node : Item {
      [Storable]
      public int index;
      [Storable]
      public double threshold;
      [Storable]
      public Node left;
      [Storable]
      public Node right;

      #region HLConstructors & Cloning
      [StorableConstructor]
      private Node(bool deserializing) : base(deserializing) { }
      private Node(Node original, Cloner cloner) : base(original, cloner) {
        index = original.index;
        threshold = original.threshold;
        left = (Node)original.left.Clone(cloner);
        right = (Node)original.right.Clone(cloner);
      }
      internal Node() {
        index = 0;
        threshold = 0;
        left = null;
        right = null;
      }
      public override IDeepCloneable Clone(Cloner cloner) {
        return new Node(this, cloner);
      }
      #endregion
    }

    [StorableClass]
    private sealed class HeapItem : Item, IComparable<HeapItem> {
      [Storable]
      public int Index;
      [Storable]
      public double Dist;

      #region HLConstructors & Cloning
      [StorableConstructor]
      private HeapItem(bool deserializing) : base(deserializing) { }
      private HeapItem(HeapItem original, Cloner cloner)
      : base(original, cloner) {
        Index = original.Index;
        Dist = original.Dist;
      }
      public override IDeepCloneable Clone(Cloner cloner) { return new HeapItem(this, cloner); }
      public HeapItem(int index, double dist) {
        Index = index;
        Dist = dist;
      }
      #endregion

      public int CompareTo(HeapItem other) {
        return Dist.CompareTo(Dist);
      }
    }
  }
}