Index: branches/RBFRegression/HeuristicLab.Algorithms.DataAnalysis/3.4/CrossValidation.cs
===================================================================
--- branches/RBFRegression/HeuristicLab.Algorithms.DataAnalysis/3.4/CrossValidation.cs	(revision 14869)
+++ branches/RBFRegression/HeuristicLab.Algorithms.DataAnalysis/3.4/CrossValidation.cs	(revision 14870)
@@ -33,4 +33,5 @@
 using HeuristicLab.Problems.DataAnalysis;
 using HeuristicLab.Problems.DataAnalysis.Symbolic;
+using HeuristicLab.Random;
 
 namespace HeuristicLab.Algorithms.DataAnalysis {
@@ -56,4 +57,5 @@
       samplesStart = new IntValue(0);
       samplesEnd = new IntValue(0);
+      shuffleSamples = new BoolValue(false);
       storeAlgorithmInEachRun = false;
 
@@ -89,4 +91,5 @@
       samplesStart = cloner.Clone(original.samplesStart);
       samplesEnd = cloner.Clone(original.samplesEnd);
+      shuffleSamples = cloner.Clone(original.shuffleSamples);
       RegisterEvents();
       if (Algorithm != null) RegisterAlgorithmEvents();
@@ -170,5 +173,9 @@
       get { return results; }
     }
-
+    [Storable]
+    private BoolValue shuffleSamples;
+    public BoolValue ShuffleSamples {
+      get { return shuffleSamples; }
+    }
     [Storable]
     private IntValue folds;
@@ -274,7 +281,14 @@
         if (clonedAlgorithms.Count == 0) {
           int testSamplesCount = (SamplesEnd.Value - SamplesStart.Value) / Folds.Value;
-
+          IDataset shuffledDataset = null;
           for (int i = 0; i < Folds.Value; i++) {
-            IAlgorithm clonedAlgorithm = (IAlgorithm)algorithm.Clone();
+            var cloner = new Cloner();
+            if (ShuffleSamples.Value) {
+              var dataAnalysisProblem = (IDataAnalysisProblem)algorithm.Problem;
+              var dataset = (Dataset)dataAnalysisProblem.ProblemData.Dataset;
+              shuffledDataset = shuffledDataset ?? dataset.Shuffle(new FastRandom());
+              cloner.RegisterClonedObject(dataset, shuffledDataset);
+            }
+            IAlgorithm clonedAlgorithm = cloner.Clone(Algorithm);
             clonedAlgorithm.Name = algorithm.Name + " Fold " + i;
             IDataAnalysisProblem problem = clonedAlgorithm.Problem as IDataAnalysisProblem;
@@ -590,9 +604,11 @@
     private void Algorithm_ExecutionStateChanged(object sender, EventArgs e) {
       switch (Algorithm.ExecutionState) {
-        case ExecutionState.Prepared: OnPrepared();
+        case ExecutionState.Prepared:
+          OnPrepared();
           break;
         case ExecutionState.Started: throw new InvalidOperationException("Algorithm template can not be started.");
         case ExecutionState.Paused: throw new InvalidOperationException("Algorithm template can not be paused.");
-        case ExecutionState.Stopped: OnStopped();
+        case ExecutionState.Stopped:
+          OnStopped();
           break;
       }
@@ -724,4 +740,5 @@
       AggregateResultValues(collectedResults);
       results.AddRange(collectedResults.Select(x => new Result(x.Key, x.Value)).Cast<IResult>().ToArray());
+      clonedAlgorithms.Clear();
       runsCounter++;
       runs.Add(new Run(string.Format("{0} Run {1}", Name, runsCounter), this));
Index: branches/RBFRegression/HeuristicLab.Algorithms.DataAnalysis/3.4/HeuristicLab.Algorithms.DataAnalysis-3.4.csproj
===================================================================
--- branches/RBFRegression/HeuristicLab.Algorithms.DataAnalysis/3.4/HeuristicLab.Algorithms.DataAnalysis-3.4.csproj	(revision 14869)
+++ branches/RBFRegression/HeuristicLab.Algorithms.DataAnalysis/3.4/HeuristicLab.Algorithms.DataAnalysis-3.4.csproj	(revision 14870)
@@ -406,4 +406,17 @@
     </Compile>
     <Compile Include="TimeSeries\AutoregressiveModeling.cs" />
+    <Compile Include="TSNE\Distances\DistanceBase.cs" />
+    <Compile Include="TSNE\Distances\EuclideanDistance.cs" />
+    <Compile Include="TSNE\Distances\IndexedItemDistance.cs" />
+    <Compile Include="TSNE\Distances\InnerProductDistance.cs" />
+    <Compile Include="TSNE\Interfaces\IDistance.cs" />
+    <Compile Include="TSNE\Interfaces\ISpacePartitioningTree.cs" />
+    <Compile Include="TSNE\Interfaces\IVantagePointTree.cs" />
+    <Compile Include="TSNE\PriorityQueue.cs" />
+    <Compile Include="TSNE\SpacePartitioningTree.cs" />
+    <Compile Include="TSNE\TSNEAlgorithm.cs" />
+    <Compile Include="TSNE\TSNEStatic.cs" />
+    <Compile Include="TSNE\TSNEUtils.cs" />
+    <Compile Include="TSNE\VantagePointTree.cs" />
   </ItemGroup>
   <ItemGroup>
Index: branches/RBFRegression/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/Distances/DistanceBase.cs
===================================================================
--- branches/RBFRegression/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/Distances/DistanceBase.cs	(revision 14870)
+++ branches/RBFRegression/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/Distances/DistanceBase.cs	(revision 14870)
@@ -0,0 +1,58 @@
+﻿#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/>.
+ */
+#endregion
+
+using System.Collections.Generic;
+using HeuristicLab.Common;
+using HeuristicLab.Core;
+using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
+
+namespace HeuristicLab.Algorithms.DataAnalysis {
+  [StorableClass]
+  public abstract class DistanceBase<T> : Item, IDistance<T> {
+
+    #region HLConstructors & Boilerplate
+    [StorableConstructor]
+    protected DistanceBase(bool deserializing) : base(deserializing) { }
+    protected DistanceBase(DistanceBase<T> original, Cloner cloner) : base(original, cloner) { }
+    protected DistanceBase() { }
+    #endregion
+
+    public abstract double Get(T a, T b);
+
+    public IComparer<T> GetDistanceComparer(T item) {
+      return new DistanceComparer(item, this);
+    }
+
+    private class DistanceComparer : IComparer<T> {
+      private readonly T item;
+      private readonly IDistance<T> dist;
+
+      public DistanceComparer(T item, IDistance<T> dist) {
+        this.dist = dist;
+        this.item = item;
+      }
+
+      public int Compare(T x, T y) {
+        return dist.Get(x, item).CompareTo(dist.Get(y, item));
+      }
+    }
+  }
+}
Index: branches/RBFRegression/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/Distances/EuclideanDistance.cs
===================================================================
--- branches/RBFRegression/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/Distances/EuclideanDistance.cs	(revision 14870)
+++ branches/RBFRegression/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/Distances/EuclideanDistance.cs	(revision 14870)
@@ -0,0 +1,51 @@
+﻿#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/>.
+ */
+#endregion
+
+using System;
+using System.Collections.Generic;
+using System.Linq;
+using HeuristicLab.Common;
+using HeuristicLab.Core;
+using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
+
+namespace HeuristicLab.Algorithms.DataAnalysis {
+  [StorableClass]
+  [Item("EuclideanDistance", "A norm function that uses Euclidean distance")]
+  public class EuclideanDistance : DistanceBase<IEnumerable<double>> {
+
+    #region HLConstructors & Boilerplate
+    [StorableConstructor]
+    protected EuclideanDistance(bool deserializing) : base(deserializing) { }
+    protected EuclideanDistance(EuclideanDistance original, Cloner cloner) : base(original, cloner) { }
+    public override IDeepCloneable Clone(Cloner cloner) { return new EuclideanDistance(this, cloner); }
+    public EuclideanDistance() { }
+    #endregion
+
+    public static double GetDistance(double[] point1, double[] point2) {
+      if (point1.Length != point2.Length) throw new ArgumentException("Euclidean distance not defined on vectors of different length");
+      return Math.Sqrt(point1.Zip(point2, (a1, b1) => (a1 - b1) * (a1 - b1)).Sum());
+    }
+
+    public override double Get(IEnumerable<double> a, IEnumerable<double> b) {
+      return GetDistance(a.ToArray(), b.ToArray());
+    }
+  }
+}
Index: branches/RBFRegression/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/Distances/IndexedItemDistance.cs
===================================================================
--- branches/RBFRegression/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/Distances/IndexedItemDistance.cs	(revision 14870)
+++ branches/RBFRegression/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/Distances/IndexedItemDistance.cs	(revision 14870)
@@ -0,0 +1,51 @@
+﻿#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/>.
+ */
+#endregion
+
+using HeuristicLab.Collections;
+using HeuristicLab.Common;
+using HeuristicLab.Core;
+using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
+
+namespace HeuristicLab.Algorithms.DataAnalysis {
+  [StorableClass]
+  [Item("IndexedItemDistance", "A distance wrapper for indexed items")]
+  internal class IndexedItemDistance<T> : DistanceBase<IndexedItem<T>> {
+    [Storable]
+    private readonly IDistance<T> dist;
+
+    #region HLConstructors & Cloning
+    [StorableConstructor]
+    protected IndexedItemDistance(bool deserializing) : base(deserializing) { }
+    protected IndexedItemDistance(IndexedItemDistance<T> original, Cloner cloner) : base(original, cloner) {
+      dist = cloner.Clone(original.dist);
+    }
+    public override IDeepCloneable Clone(Cloner cloner) { return new IndexedItemDistance<T>(this, cloner); }
+    public IndexedItemDistance(IDistance<T> distance) {
+      dist = distance;
+    }
+    #endregion
+
+    public override double Get(IndexedItem<T> a, IndexedItem<T> b) {
+      return dist.Get(a.Value, b.Value);
+    }
+  }
+}
+
Index: branches/RBFRegression/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/Distances/InnerProductDistance.cs
===================================================================
--- branches/RBFRegression/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/Distances/InnerProductDistance.cs	(revision 14870)
+++ branches/RBFRegression/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/Distances/InnerProductDistance.cs	(revision 14870)
@@ -0,0 +1,67 @@
+﻿#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/>.
+ */
+#endregion
+
+using System;
+using System.Collections.Generic;
+using System.Linq;
+using HeuristicLab.Common;
+using HeuristicLab.Core;
+using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
+
+namespace HeuristicLab.Algorithms.DataAnalysis {
+
+  /// <summary>
+  /// The angluar distance as defined as a normalized distance measure dependent on the angle between two vectors.
+  /// It is designed for vectors with all positive coordinates.
+  /// </summary>
+  [StorableClass]
+  [Item("InnerProductDistance", "The angluar distance as defined as a normalized distance measure dependent on the angle between two vectors.\nIt is designed for vectors with all positive coordinates")]
+  public class InnerProductDistance : DistanceBase<IEnumerable<double>> {
+
+    #region HLConstructors
+    [StorableConstructor]
+    protected InnerProductDistance(bool deserializing) : base(deserializing) { }
+    protected InnerProductDistance(InnerProductDistance original, Cloner cloner)
+      : base(original, cloner) { }
+    public InnerProductDistance() { }
+    public override IDeepCloneable Clone(Cloner cloner) {
+      return new InnerProductDistance(this, cloner);
+    }
+    #endregion
+
+    #region statics
+    public static double GetDistance(IEnumerable<double> point1, IEnumerable<double> point2) {
+      var xs = point1.GetEnumerator();
+      var ys = point2.GetEnumerator();
+      var sum = 0.0;
+      while(xs.MoveNext() & ys.MoveNext()) {
+        if(xs.Current < 0 || ys.Current < 0) throw new ArgumentException("Inner product distance is only defined for vectors with non-negative elements");
+        sum += xs.Current * ys.Current;
+      }
+      if(xs.MoveNext() || ys.MoveNext()) throw new ArgumentException("Enumerables contain a different number of elements");
+      return sum;
+    }
+    #endregion
+    public override double Get(IEnumerable<double> a, IEnumerable<double> b) {
+      return GetDistance(a, b);
+    }
+  }
+}
Index: branches/RBFRegression/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/Interfaces/IDistance.cs
===================================================================
--- branches/RBFRegression/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/Interfaces/IDistance.cs	(revision 14870)
+++ branches/RBFRegression/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/Interfaces/IDistance.cs	(revision 14870)
@@ -0,0 +1,43 @@
+﻿#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/>.
+ */
+#endregion
+
+using System.Collections.Generic;
+using HeuristicLab.Core;
+
+namespace HeuristicLab.Algorithms.DataAnalysis {
+  public interface IDistance<in T> : IItem {
+    /// <summary>
+    /// Calculates a distance measure between two objects.
+    /// 1.) non-negative d(x,y) >= 0
+    /// 2.) symmetric d(x,y) = d(y,x)
+    /// 3.) zero-reflexive d(x,x) = 0;
+    /// </summary>
+    /// <returns>d(x,y)</returns>
+    double Get(T x, T y);
+
+    /// <summary>
+    /// Returns a comparator wich compares the distances to item. (allows for sorting nearest/farthest neighbours)
+    /// </summary>
+    /// <param name="item"></param>
+    /// <returns></returns>
+    IComparer<T> GetDistanceComparer(T item);
+  }
+}
Index: branches/RBFRegression/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/Interfaces/ISpacePartitioningTree.cs
===================================================================
--- branches/RBFRegression/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/Interfaces/ISpacePartitioningTree.cs	(revision 14870)
+++ branches/RBFRegression/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/Interfaces/ISpacePartitioningTree.cs	(revision 14870)
@@ -0,0 +1,37 @@
+﻿#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/>.
+ */
+#endregion
+
+using HeuristicLab.Common;
+
+namespace HeuristicLab.Algorithms.DataAnalysis {
+  public interface ISpacePartitioningTree {
+
+    int GetDepth();
+    bool Insert(int newIndex);
+    void Subdivide();
+    bool IsCorrect();
+    void GetAllIndices(int[] indices);
+    int GetAllIndices(int[] indices, int loc);
+    void ComputeNonEdgeForces(int pointIndex, double theta, double[] negF, ref double sumQ);
+    void ComputeEdgeForces(int[] rowP, int[] colP, double[] valP, int n, double[,] posF);
+
+  }
+}
Index: branches/RBFRegression/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/Interfaces/IVantagePointTree.cs
===================================================================
--- branches/RBFRegression/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/Interfaces/IVantagePointTree.cs	(revision 14870)
+++ branches/RBFRegression/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/Interfaces/IVantagePointTree.cs	(revision 14870)
@@ -0,0 +1,29 @@
+﻿#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/>.
+ */
+#endregion
+
+using System.Collections.Generic;
+using HeuristicLab.Common;
+
+namespace HeuristicLab.Algorithms.DataAnalysis {
+  public interface IVantagePointTree<T> {
+    void Search(T target, int k, out IList<T> results, out IList<double> distances);
+  }
+}
Index: branches/RBFRegression/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/PriorityQueue.cs
===================================================================
--- branches/RBFRegression/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/PriorityQueue.cs	(revision 14870)
+++ branches/RBFRegression/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/PriorityQueue.cs	(revision 14870)
@@ -0,0 +1,202 @@
+﻿#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/>.
+ */
+#endregion
+
+using System;
+using System.Collections.Generic;
+using System.Linq;
+
+namespace HeuristicLab.Algorithms.DataAnalysis {
+  // TODO: move to HeuristicLab.Collections
+
+  // Code Taken from branch: RoutePlanning, Heap4
+  // implementation based on C++ version from Peter Sanders
+  // http://www.mpi-inf.mpg.de/~sanders/programs/spq/
+  public sealed class PriorityQueue<TK, TV> where TK : IComparable {
+    private class KNElement {
+      public TK Key { get; set; }
+      public TV Value { get; set; }
+    }
+
+    private readonly int capacity;
+    private int size; // index of last used element
+    private int finalLayerSize; // size of first layer with free space
+    private int finalLayerDist; // distance to end of layer
+                                //private KNElement[] rawData;
+    private readonly KNElement[] data;    // aligned version of rawData
+
+    public PriorityQueue(TK supremum, TK infimum, int cap) {
+      capacity = cap;
+      data = Enumerable.Range(0, capacity + 4).Select(x => new KNElement()).ToArray();
+
+      data[0].Key = infimum; // sentinel
+      data[capacity + 1].Key = supremum;
+      data[capacity + 2].Key = supremum;
+      data[capacity + 3].Key = supremum;
+      Reset();
+    }
+
+    public int Size {
+      get { return size; }
+    }
+
+    public TK Supremum {
+      get { return data[capacity + 1].Key; }
+    }
+
+    public KeyValuePair<TK, TV> PeekMin() {
+      if(size == 0) {
+        throw new InvalidOperationException("Heap is empty");
+      }
+      return new KeyValuePair<TK, TV>(data[1].Key, data[1].Value);
+    }
+
+    public TK PeekMinKey() {
+      return data[1].Key;
+    }
+
+    public TV PeekMinValue() {
+      return data[1].Value;
+    }
+
+    public void RemoveMin() {
+      // first move up elements on a min-path
+      var hole = 1;
+      var succ = 2;
+      var layerSize = 4;
+      var layerPos = 0;
+      var sz = size;
+      size = sz - 1;
+      finalLayerDist++;
+      if(finalLayerDist == finalLayerSize) {
+        finalLayerSize >>= 2;
+        finalLayerDist = 0;
+      }
+
+      while(succ < sz) {
+        var minKey = data[succ].Key;
+        var delta = 0;
+
+        var otherKey = data[succ + 1].Key;
+        if(otherKey.CompareTo(minKey) < 0) {
+          minKey = otherKey;
+          delta = 1;
+        }
+        otherKey = data[succ + 2].Key;
+        if(otherKey.CompareTo(minKey) < 0) {
+          minKey = otherKey;
+          delta = 2;
+        }
+        otherKey = data[succ + 3].Key;
+        if(otherKey.CompareTo(minKey) < 0) {
+          minKey = otherKey;
+          delta = 3;
+        }
+        succ += delta;
+        layerPos += delta;
+
+        // move min successor up
+        data[hole].Key = minKey;
+        data[hole].Value = data[succ].Value;
+
+        // step to next layer
+        hole = succ;
+        succ = succ - layerPos + layerSize; // beginning of next layer
+        layerPos <<= 2;
+        succ += layerPos; // now correct value
+        layerSize <<= 2;
+      }
+
+      // bubble up rightmost element
+      var bubble = data[sz].Key;
+      layerSize >>= 2; // now size of hole's layer
+      layerPos >>= 2; // now pos of hole within its layer
+
+      var layerDist = layerSize - layerPos - 1; // hole's dist to end of layer
+      var pred = hole + layerDist - layerSize; // end of pred's layer for now
+      layerSize >>= 2; // now size of pred's layer
+      layerDist >>= 2; // now pred's pos in layer
+      pred = pred - layerDist; // finally preds index
+
+      while(data[pred].Key.CompareTo(bubble) > 0) {  // must terminate since inf at root
+        data[hole].Key = data[pred].Key;
+        data[hole].Value = data[pred].Value;
+        hole = pred;
+        pred = hole + layerDist - layerSize; // end of hole's layer for now
+        layerSize >>= 2;
+        layerDist >>= 2;
+        pred = pred - layerDist; // finally preds index
+      }
+
+      // finally move data to hole
+      data[hole].Key = bubble;
+      data[hole].Value = data[sz].Value;
+
+      data[sz].Key = Supremum; // mark as deleted
+    }
+
+    public void Insert(TK key, TV value) {
+      var layerSize = finalLayerSize;
+      var layerDist = finalLayerDist;
+      finalLayerDist--;
+
+      if(finalLayerDist == -1) { // layer full
+                                 // start next layer
+        finalLayerSize <<= 2;
+        finalLayerDist = finalLayerSize - 1;
+      }
+
+      size++;
+
+      var hole = size;
+      var pred = hole + layerDist - layerSize; // end of pred's layer for now
+      layerSize >>= 2; // now size of pred's layer
+      layerDist >>= 2;
+      pred = pred - layerDist; // finally preds index
+      var predKey = data[pred].Key;
+      while(predKey.CompareTo(key) > 0) {
+        data[hole].Key = predKey;
+        data[hole].Value = data[pred].Value;
+        hole = pred;
+        pred = hole + layerDist - layerSize; // end of pred's layer for now
+        layerSize >>= 2; // now isze of pred's layer
+        layerDist >>= 2;
+        pred = pred - layerDist; // finally preds index
+        predKey = data[pred].Key;
+      }
+
+      // finally move data to hole
+      data[hole].Key = key;
+      data[hole].Value = value;
+    }
+
+    // reset size to 0 and fill data array with sentinels
+    private void Reset() {
+      size = 0;
+      finalLayerSize = 1;
+      finalLayerDist = 0;
+      var sup = Supremum;
+      var cap = capacity;
+      for(var i = 1; i <= cap; i++) {
+        data[i].Key = sup;
+      }
+    }
+  }
+}
Index: branches/RBFRegression/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/SpacePartitioningTree.cs
===================================================================
--- branches/RBFRegression/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/SpacePartitioningTree.cs	(revision 14870)
+++ branches/RBFRegression/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/SpacePartitioningTree.cs	(revision 14870)
@@ -0,0 +1,326 @@
+﻿#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;
+
+namespace HeuristicLab.Algorithms.DataAnalysis {
+  /// <summary>
+  /// Space partitioning tree (SPTree)
+  /// </summary>
+  public class SpacePartitioningTree : ISpacePartitioningTree {
+    private const uint QT_NODE_CAPACITY = 1;
+
+    private double[] buff;
+    private SpacePartitioningTree parent;
+    private int dimension;
+    private bool isLeaf;
+    private uint size;
+    private uint cumulativeSize;
+
+    // Axis-aligned bounding box stored as a center with half-dimensions to represent the boundaries of this quad tree
+    private Cell boundary;
+
+    private double[,] data;
+
+    // Indices in this space-partitioning tree node, corresponding center-of-mass, and list of all children
+    private double[] centerOfMass;
+    private readonly int[] index = new int[QT_NODE_CAPACITY];
+
+    // Children
+    private SpacePartitioningTree[] children;
+    private uint noChildren;
+
+    public SpacePartitioningTree(double[,] inpData) {
+      var d = inpData.GetLength(1);
+      var n = inpData.GetLength(0);
+      var meanY = new double[d];
+      var minY = new double[d];
+      for (var i = 0; i < d; i++) minY[i] = double.MaxValue;
+      var maxY = new double[d];
+      for (var i = 0; i < d; i++) maxY[i] = double.MinValue;
+      for (uint i = 0; i < n; i++) {
+        for (uint j = 0; j < d; j++) {
+          meanY[j] += inpData[i, j];
+          if (inpData[i, j] < minY[j]) minY[j] = inpData[i, j];
+          if (inpData[i, j] > maxY[j]) maxY[j] = inpData[i, j];
+        }
+      }
+      for (var i = 0; i < d; i++) meanY[i] /= n;
+      var width = new double[d];
+      for (var i = 0; i < d; i++) width[i] = Math.Max(maxY[i] - meanY[i], meanY[i] - minY[i]) + 1e-5;
+      Init(null, inpData, meanY, width);
+      Fill(n);
+    }
+
+    public SpacePartitioningTree(double[,] inpData, IEnumerable<double> impCorner, IEnumerable<double> impWith) {
+      Init(null, inpData, impCorner, impWith);
+    }
+    public SpacePartitioningTree(SpacePartitioningTree parent, double[,] inpData, IEnumerable<double> impCorner, IEnumerable<double> impWith) {
+      Init(parent, inpData, impCorner, impWith);
+    }
+
+    public ISpacePartitioningTree GetParent() {
+      return parent;
+    }
+
+    public bool Insert(int newIndex) {
+      // Ignore objects which do not belong in this quad tree
+      var point = new double[dimension];
+      Buffer.BlockCopy(data, sizeof(double) * dimension * newIndex, point, 0, sizeof(double) * dimension);
+      if (!boundary.ContainsPoint(point)) return false;
+      cumulativeSize++;
+      // Online update of cumulative size and center-of-mass
+      var mult1 = (double)(cumulativeSize - 1) / cumulativeSize;
+      var mult2 = 1.0 / cumulativeSize;
+      for (var i = 0; i < dimension; i++) centerOfMass[i] *= mult1;
+      for (var i = 0; i < dimension; i++) centerOfMass[i] += mult2 * point[i];
+
+      // If there is space in this quad tree and it is a leaf, add the object here
+      if (isLeaf && size < QT_NODE_CAPACITY) {
+        index[size] = newIndex;
+        size++;
+        return true;
+      }
+
+      // Don't add duplicates for now (this is not very nice)
+      var anyDuplicate = false;
+      for (uint n = 0; n < size; n++) {
+        var duplicate = true;
+        for (var d = 0; d < dimension; d++) {
+          if (Math.Abs(point[d] - data[index[n], d]) < double.Epsilon) continue;
+          duplicate = false; break;
+        }
+        anyDuplicate = anyDuplicate | duplicate;
+      }
+      if (anyDuplicate) return true;
+
+      // Otherwise, we need to subdivide the current cell
+      if (isLeaf) Subdivide();
+      // Find out where the point can be inserted
+      for (var i = 0; i < noChildren; i++) {
+        if (children[i].Insert(newIndex)) return true;
+      }
+
+      // Otherwise, the point cannot be inserted (this should never happen)
+      return false;
+    }
+
+    public void Subdivide() {
+      // Create new children
+      var newCorner = new double[dimension];
+      var newWidth = new double[dimension];
+      for (var i = 0; i < noChildren; i++) {
+        var div = 1;
+        for (var d = 0; d < dimension; d++) {
+          newWidth[d] = .5 * boundary.GetWidth(d);
+          if ((i / div) % 2 == 1) newCorner[d] = boundary.GetCorner(d) - .5 * boundary.GetWidth(d);
+          else newCorner[d] = boundary.GetCorner(d) + .5 * boundary.GetWidth(d);
+          div *= 2;
+        }
+        children[i] = new SpacePartitioningTree(this, data, newCorner, newWidth);
+      }
+
+      // Move existing points to correct children
+      for (var i = 0; i < size; i++) {
+        var success = false;
+        for (var j = 0; j < noChildren; j++) {
+          if (!success) success = children[j].Insert(index[i]);
+        }
+        index[i] = -1; // as in tSNE implementation by van der Maaten
+      }
+
+      // Empty parent node
+      size = 0;
+      isLeaf = false;
+    }
+
+    public bool IsCorrect() {
+      var row = new double[dimension];
+      for (var n = 0; n < size; n++) {
+        Buffer.BlockCopy(data, sizeof(double) * dimension * index[n], row, 0, sizeof(double) * dimension);
+        if (!boundary.ContainsPoint(row)) return false;
+      }
+      if (isLeaf) return true;
+      var correct = true;
+      for (var i = 0; i < noChildren; i++) correct = correct && children[i].IsCorrect();
+      return correct;
+    }
+
+    public void GetAllIndices(int[] indices) {
+      GetAllIndices(indices, 0);
+    }
+
+    public int GetAllIndices(int[] indices, int loc) {
+      // Gather indices in current quadrant
+      for (var i = 0; i < size; i++) indices[loc + i] = index[i];
+      loc += (int)size;
+      // Gather indices in children
+      if (isLeaf) return loc;
+      for (var i = 0; i < noChildren; i++) loc = children[i].GetAllIndices(indices, loc);
+      return loc;
+    }
+
+    public int GetDepth() {
+      return isLeaf ? 1 : 1 + children.Max(x => x.GetDepth());
+    }
+
+    public void ComputeNonEdgeForces(int pointIndex, double theta, double[] negF, ref double sumQ) {
+      // Make sure that we spend no time on empty nodes or self-interactions
+      if (cumulativeSize == 0 || (isLeaf && size == 1 && index[0] == pointIndex)) return;
+
+      // Compute distance between point and center-of-mass
+      var D = .0;
+      for (var d = 0; d < dimension; d++) buff[d] = data[pointIndex, d] - centerOfMass[d];
+      for (var d = 0; d < dimension; d++) D += buff[d] * buff[d];
+
+      // Check whether we can use this node as a "summary"
+      var maxWidth = 0.0;
+      for (var d = 0; d < dimension; d++) {
+        var curWidth = boundary.GetWidth(d);
+        maxWidth = (maxWidth > curWidth) ? maxWidth : curWidth;
+      }
+      if (isLeaf || maxWidth / Math.Sqrt(D) < theta) {
+
+        // Compute and add t-SNE force between point and current node
+        D = 1.0 / (1.0 + D);
+        var mult = cumulativeSize * D;
+        sumQ += mult;
+        mult *= D;
+        for (var d = 0; d < dimension; d++) negF[d] += mult * buff[d];
+      } else {
+
+        // Recursively apply Barnes-Hut to children
+        for (var i = 0; i < noChildren; i++) children[i].ComputeNonEdgeForces(pointIndex, theta, negF, ref sumQ);
+      }
+    }
+
+    // does not use the tree
+    public void ComputeEdgeForces(int[] rowP, int[] colP, double[] valP, int n, double[,] posF) {
+      // Loop over all edges in the graph
+      for (var k = 0; k < n; k++) {
+        for (var i = rowP[k]; i < rowP[k + 1]; i++) {
+
+          // Compute pairwise distance and Q-value
+          // uses squared distance
+          var d = 1.0;
+          for (var j = 0; j < dimension; j++) buff[j] = data[k, j] - data[colP[i], j];
+          for (var j = 0; j < dimension; j++) d += buff[j] * buff[j];
+          d = valP[i] / d;
+
+          // Sum positive force
+          for (var j = 0; j < dimension; j++) posF[k, j] += d * buff[j];
+        }
+      }
+    }
+
+    #region Helpers
+    private void Fill(int n) {
+      for (var i = 0; i < n; i++) Insert(i);
+    }
+    private void Init(SpacePartitioningTree p, double[,] inpData, IEnumerable<double> inpCorner, IEnumerable<double> inpWidth) {
+      parent = p;
+      dimension = inpData.GetLength(1);
+      noChildren = 2;
+      for (uint i = 1; i < dimension; i++) noChildren *= 2;
+      data = inpData;
+      isLeaf = true;
+      size = 0;
+      cumulativeSize = 0;
+      boundary = new Cell((uint)dimension);
+      inpCorner.ForEach((i, x) => boundary.SetCorner(i, x));
+      inpWidth.ForEach((i, x) => boundary.SetWidth(i, x));
+
+      children = new SpacePartitioningTree[noChildren];
+      centerOfMass = new double[dimension];
+      buff = new double[dimension];
+
+    }
+    #endregion
+
+
+    private class Cell {
+      private readonly uint dimension;
+      private readonly double[] corner;
+      private readonly double[] width;
+
+      public Cell(uint inpDimension) {
+        dimension = inpDimension;
+        corner = new double[dimension];
+        width = new double[dimension];
+      }
+
+      public double GetCorner(int d) {
+        return corner[d];
+      }
+      public double GetWidth(int d) {
+        return width[d];
+      }
+      public void SetCorner(int d, double val) {
+        corner[d] = val;
+      }
+      public void SetWidth(int d, double val) {
+        width[d] = val;
+      }
+      public bool ContainsPoint(double[] point) {
+        for (var d = 0; d < dimension; d++)
+          if (corner[d] - width[d] > point[d] || corner[d] + width[d] < point[d]) return false;
+        return true;
+      }
+    }
+  }
+}
Index: branches/RBFRegression/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/TSNEAlgorithm.cs
===================================================================
--- branches/RBFRegression/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/TSNEAlgorithm.cs	(revision 14870)
+++ branches/RBFRegression/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/TSNEAlgorithm.cs	(revision 14870)
@@ -0,0 +1,419 @@
+#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/>.
+ */
+#endregion
+
+using System;
+using System.Collections.Generic;
+using System.Drawing;
+using System.Linq;
+using System.Threading;
+using HeuristicLab.Analysis;
+using HeuristicLab.Common;
+using HeuristicLab.Core;
+using HeuristicLab.Data;
+using HeuristicLab.Optimization;
+using HeuristicLab.Parameters;
+using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
+using HeuristicLab.Problems.DataAnalysis;
+using HeuristicLab.Random;
+
+namespace HeuristicLab.Algorithms.DataAnalysis {
+  /// <summary>
+  /// t-distributed stochastic neighbourhood embedding (tSNE) projects the data in a low dimensional 
+  /// space to allow visual cluster identification.
+  /// </summary>
+  [Item("tSNE", "t-distributed stochastic neighbourhood embedding projects the data in a low " +
+                "dimensional space to allow visual cluster identification.")]
+  [Creatable(CreatableAttribute.Categories.DataAnalysis, Priority = 100)]
+  [StorableClass]
+  public sealed class TSNEAlgorithm : BasicAlgorithm {
+    public override bool SupportsPause {
+      get { return true; }
+    }
+    public override Type ProblemType {
+      get { return typeof(IDataAnalysisProblem); }
+    }
+    public new IDataAnalysisProblem Problem {
+      get { return (IDataAnalysisProblem)base.Problem; }
+      set { base.Problem = value; }
+    }
+
+    #region parameter names
+    private const string DistanceParameterName = "DistanceFunction";
+    private const string PerplexityParameterName = "Perplexity";
+    private const string ThetaParameterName = "Theta";
+    private const string NewDimensionsParameterName = "Dimensions";
+    private const string MaxIterationsParameterName = "MaxIterations";
+    private const string StopLyingIterationParameterName = "StopLyingIteration";
+    private const string MomentumSwitchIterationParameterName = "MomentumSwitchIteration";
+    private const string InitialMomentumParameterName = "InitialMomentum";
+    private const string FinalMomentumParameterName = "FinalMomentum";
+    private const string EtaParameterName = "Eta";
+    private const string SetSeedRandomlyParameterName = "SetSeedRandomly";
+    private const string SeedParameterName = "Seed";
+    private const string ClassesParameterName = "ClassNames";
+    private const string NormalizationParameterName = "Normalization";
+    private const string UpdateIntervalParameterName = "UpdateInterval";
+    #endregion
+
+    #region result names
+    private const string IterationResultName = "Iteration";
+    private const string ErrorResultName = "Error";
+    private const string ErrorPlotResultName = "Error plot";
+    private const string ScatterPlotResultName = "Scatterplot";
+    private const string DataResultName = "Projected data";
+    #endregion
+
+    #region parameter properties
+    public IFixedValueParameter<DoubleValue> PerplexityParameter {
+      get { return Parameters[PerplexityParameterName] as IFixedValueParameter<DoubleValue>; }
+    }
+    public IFixedValueParameter<DoubleValue> ThetaParameter {
+      get { return Parameters[ThetaParameterName] as IFixedValueParameter<DoubleValue>; }
+    }
+    public IFixedValueParameter<IntValue> NewDimensionsParameter {
+      get { return Parameters[NewDimensionsParameterName] as IFixedValueParameter<IntValue>; }
+    }
+    public IValueParameter<IDistance<double[]>> DistanceParameter {
+      get { return Parameters[DistanceParameterName] as IValueParameter<IDistance<double[]>>; }
+    }
+    public IFixedValueParameter<IntValue> MaxIterationsParameter {
+      get { return Parameters[MaxIterationsParameterName] as IFixedValueParameter<IntValue>; }
+    }
+    public IFixedValueParameter<IntValue> StopLyingIterationParameter {
+      get { return Parameters[StopLyingIterationParameterName] as IFixedValueParameter<IntValue>; }
+    }
+    public IFixedValueParameter<IntValue> MomentumSwitchIterationParameter {
+      get { return Parameters[MomentumSwitchIterationParameterName] as IFixedValueParameter<IntValue>; }
+    }
+    public IFixedValueParameter<DoubleValue> InitialMomentumParameter {
+      get { return Parameters[InitialMomentumParameterName] as IFixedValueParameter<DoubleValue>; }
+    }
+    public IFixedValueParameter<DoubleValue> FinalMomentumParameter {
+      get { return Parameters[FinalMomentumParameterName] as IFixedValueParameter<DoubleValue>; }
+    }
+    public IFixedValueParameter<DoubleValue> EtaParameter {
+      get { return Parameters[EtaParameterName] as IFixedValueParameter<DoubleValue>; }
+    }
+    public IFixedValueParameter<BoolValue> SetSeedRandomlyParameter {
+      get { return Parameters[SetSeedRandomlyParameterName] as IFixedValueParameter<BoolValue>; }
+    }
+    public IFixedValueParameter<IntValue> SeedParameter {
+      get { return Parameters[SeedParameterName] as IFixedValueParameter<IntValue>; }
+    }
+    public IFixedValueParameter<StringValue> ClassesParameter {
+      get { return Parameters[ClassesParameterName] as IFixedValueParameter<StringValue>; }
+    }
+    public IFixedValueParameter<BoolValue> NormalizationParameter {
+      get { return Parameters[NormalizationParameterName] as IFixedValueParameter<BoolValue>; }
+    }
+    public IFixedValueParameter<IntValue> UpdateIntervalParameter {
+      get { return Parameters[UpdateIntervalParameterName] as IFixedValueParameter<IntValue>; }
+    }
+    #endregion
+
+    #region  Properties
+    public IDistance<double[]> Distance {
+      get { return DistanceParameter.Value; }
+    }
+    public double Perplexity {
+      get { return PerplexityParameter.Value.Value; }
+      set { PerplexityParameter.Value.Value = value; }
+    }
+    public double Theta {
+      get { return ThetaParameter.Value.Value; }
+      set { ThetaParameter.Value.Value = value; }
+    }
+    public int NewDimensions {
+      get { return NewDimensionsParameter.Value.Value; }
+      set { NewDimensionsParameter.Value.Value = value; }
+    }
+    public int MaxIterations {
+      get { return MaxIterationsParameter.Value.Value; }
+      set { MaxIterationsParameter.Value.Value = value; }
+    }
+    public int StopLyingIteration {
+      get { return StopLyingIterationParameter.Value.Value; }
+      set { StopLyingIterationParameter.Value.Value = value; }
+    }
+    public int MomentumSwitchIteration {
+      get { return MomentumSwitchIterationParameter.Value.Value; }
+      set { MomentumSwitchIterationParameter.Value.Value = value; }
+    }
+    public double InitialMomentum {
+      get { return InitialMomentumParameter.Value.Value; }
+      set { InitialMomentumParameter.Value.Value = value; }
+    }
+    public double FinalMomentum {
+      get { return FinalMomentumParameter.Value.Value; }
+      set { FinalMomentumParameter.Value.Value = value; }
+    }
+    public double Eta {
+      get { return EtaParameter.Value.Value; }
+      set { EtaParameter.Value.Value = value; }
+    }
+    public bool SetSeedRandomly {
+      get { return SetSeedRandomlyParameter.Value.Value; }
+      set { SetSeedRandomlyParameter.Value.Value = value; }
+    }
+    public int Seed {
+      get { return SeedParameter.Value.Value; }
+      set { SeedParameter.Value.Value = value; }
+    }
+    public string Classes {
+      get { return ClassesParameter.Value.Value; }
+      set { ClassesParameter.Value.Value = value; }
+    }
+    public bool Normalization {
+      get { return NormalizationParameter.Value.Value; }
+      set { NormalizationParameter.Value.Value = value; }
+    }
+
+    public int UpdateInterval {
+      get { return UpdateIntervalParameter.Value.Value; }
+      set { UpdateIntervalParameter.Value.Value = value; }
+    }
+    #endregion
+
+    #region Constructors & Cloning
+    [StorableConstructor]
+    private TSNEAlgorithm(bool deserializing) : base(deserializing) { }
+
+    private TSNEAlgorithm(TSNEAlgorithm original, Cloner cloner) : base(original, cloner) {
+      if(original.dataRowNames!=null)
+      this.dataRowNames = new Dictionary<string, List<int>>(original.dataRowNames);
+      if (original.dataRows != null)
+        this.dataRows = original.dataRows.ToDictionary(kvp => kvp.Key, kvp => cloner.Clone(kvp.Value));
+      if (original.state != null)
+        this.state = cloner.Clone(original.state);
+      this.iter = original.iter;
+    }
+    public override IDeepCloneable Clone(Cloner cloner) { return new TSNEAlgorithm(this, cloner); }
+    public TSNEAlgorithm() {
+      Problem = new RegressionProblem();
+      Parameters.Add(new ValueParameter<IDistance<double[]>>(DistanceParameterName, "The distance function used to differentiate similar from non-similar points", new EuclideanDistance()));
+      Parameters.Add(new FixedValueParameter<DoubleValue>(PerplexityParameterName, "Perplexity-parameter of tSNE. Comparable to k in a k-nearest neighbour algorithm. Recommended value is floor(number of points /3) or lower", new DoubleValue(25)));
+      Parameters.Add(new FixedValueParameter<DoubleValue>(ThetaParameterName, "Value describing how much appoximated " +
+                                                                              "gradients my differ from exact gradients. Set to 0 for exact calculation and in [0,1] otherwise. " +
+                                                                              "Appropriate values for theta are between 0.1 and 0.7 (default = 0.5). CAUTION: exact calculation of " +
+                                                                              "forces requires building a non-sparse N*N matrix where N is the number of data points. This may " +
+                                                                              "exceed memory limitations. The function is designed to run on large (N > 5000) data sets. It may give" +
+                                                                              " poor performance on very small data sets(it is better to use a standard t - SNE implementation on such data).", new DoubleValue(0)));
+      Parameters.Add(new FixedValueParameter<IntValue>(NewDimensionsParameterName, "Dimensionality of projected space (usually 2 for easy visual analysis)", new IntValue(2)));
+      Parameters.Add(new FixedValueParameter<IntValue>(MaxIterationsParameterName, "Maximum number of iterations for gradient descent.", new IntValue(1000)));
+      Parameters.Add(new FixedValueParameter<IntValue>(StopLyingIterationParameterName, "Number of iterations after which p is no longer approximated.", new IntValue(0)));
+      Parameters.Add(new FixedValueParameter<IntValue>(MomentumSwitchIterationParameterName, "Number of iterations after which the momentum in the gradient descent is switched.", new IntValue(0)));
+      Parameters.Add(new FixedValueParameter<DoubleValue>(InitialMomentumParameterName, "The initial momentum in the gradient descent.", new DoubleValue(0.5)));
+      Parameters.Add(new FixedValueParameter<DoubleValue>(FinalMomentumParameterName, "The final momentum.", new DoubleValue(0.8)));
+      Parameters.Add(new FixedValueParameter<DoubleValue>(EtaParameterName, "Gradient descent learning rate.", new DoubleValue(10)));
+      Parameters.Add(new FixedValueParameter<BoolValue>(SetSeedRandomlyParameterName, "If the seed should be random.", new BoolValue(true)));
+      Parameters.Add(new FixedValueParameter<IntValue>(SeedParameterName, "The seed used if it should not be random.", new IntValue(0)));
+      Parameters.Add(new FixedValueParameter<StringValue>(ClassesParameterName, "name of the column specifying the class lables of each data point. If the label column can not be found training/test is used as labels.", new StringValue("none")));
+      Parameters.Add(new FixedValueParameter<BoolValue>(NormalizationParameterName, "Whether the data should be zero centered and have variance of 1 for each variable, so different scalings are ignored.", new BoolValue(true)));
+      Parameters.Add(new FixedValueParameter<IntValue>(UpdateIntervalParameterName, "", new IntValue(50)));
+      Parameters[UpdateIntervalParameterName].Hidden = true;
+
+      MomentumSwitchIterationParameter.Hidden = true;
+      InitialMomentumParameter.Hidden = true;
+      FinalMomentumParameter.Hidden = true;
+      StopLyingIterationParameter.Hidden = true;
+      EtaParameter.Hidden = false;
+    }
+    #endregion
+
+    [Storable]
+    private Dictionary<string, List<int>> dataRowNames;
+    [Storable]
+    private Dictionary<string, ScatterPlotDataRow> dataRows;
+    [Storable]
+    private TSNEStatic<double[]>.TSNEState state;
+    [Storable]
+    private int iter;
+
+    public override void Prepare() {
+      base.Prepare();
+      dataRowNames = null;
+      dataRows = null;
+      state = null;
+    }
+
+    protected override void Run(CancellationToken cancellationToken) {
+      var problemData = Problem.ProblemData;
+      // set up and initialized everything if necessary
+      if (state == null) {
+        if (SetSeedRandomly) Seed = new System.Random().Next();
+        var random = new MersenneTwister((uint)Seed);
+        var dataset = problemData.Dataset;
+        var allowedInputVariables = problemData.AllowedInputVariables.ToArray();
+        var data = new double[dataset.Rows][];
+        for (var row = 0; row < dataset.Rows; row++)
+          data[row] = allowedInputVariables.Select(col => dataset.GetDoubleValue(col, row)).ToArray();
+
+        if (Normalization) data = NormalizeData(data);
+
+        state = TSNEStatic<double[]>.CreateState(data, Distance, random, NewDimensions, Perplexity, Theta,
+          StopLyingIteration, MomentumSwitchIteration, InitialMomentum, FinalMomentum, Eta);
+
+        SetUpResults(data);
+        iter = 0;
+      }
+      for (; iter < MaxIterations && !cancellationToken.IsCancellationRequested; iter++) {
+        if (iter % UpdateInterval == 0)
+          Analyze(state);
+        TSNEStatic<double[]>.Iterate(state);
+      }
+      Analyze(state);
+    }
+
+    private void SetUpResults(IReadOnlyCollection<double[]> data) {
+      if (Results == null) return;
+      var results = Results;
+      dataRowNames = new Dictionary<string, List<int>>();
+      dataRows = new Dictionary<string, ScatterPlotDataRow>();
+      var problemData = Problem.ProblemData;
+
+      //color datapoints acording to classes variable (be it double or string)
+      if (problemData.Dataset.VariableNames.Contains(Classes)) {
+        if ((problemData.Dataset as Dataset).VariableHasType<string>(Classes)) {
+          var classes = problemData.Dataset.GetStringValues(Classes).ToArray();
+          for (var i = 0; i < classes.Length; i++) {
+            if (!dataRowNames.ContainsKey(classes[i])) dataRowNames.Add(classes[i], new List<int>());
+            dataRowNames[classes[i]].Add(i);
+          }
+        } else if ((problemData.Dataset as Dataset).VariableHasType<double>(Classes)) {
+          var classValues = problemData.Dataset.GetDoubleValues(Classes).ToArray();
+          var max = classValues.Max() + 0.1;
+          var min = classValues.Min() - 0.1;
+          const int contours = 8;
+          for (var i = 0; i < contours; i++) {
+            var contourname = GetContourName(i, min, max, contours);
+            dataRowNames.Add(contourname, new List<int>());
+            dataRows.Add(contourname, new ScatterPlotDataRow(contourname, "", new List<Point2D<double>>()));
+            dataRows[contourname].VisualProperties.Color = GetHeatMapColor(i, contours);
+            dataRows[contourname].VisualProperties.PointSize = i + 3;
+          }
+          for (var i = 0; i < classValues.Length; i++) {
+            dataRowNames[GetContourName(classValues[i], min, max, contours)].Add(i);
+          }
+        }
+      } else {
+        dataRowNames.Add("Training", problemData.TrainingIndices.ToList());
+        dataRowNames.Add("Test", problemData.TestIndices.ToList());
+      }
+
+      if (!results.ContainsKey(IterationResultName)) results.Add(new Result(IterationResultName, new IntValue(0)));
+      else ((IntValue)results[IterationResultName].Value).Value = 0;
+
+      if (!results.ContainsKey(ErrorResultName)) results.Add(new Result(ErrorResultName, new DoubleValue(0)));
+      else ((DoubleValue)results[ErrorResultName].Value).Value = 0;
+
+      if (!results.ContainsKey(ErrorPlotResultName)) results.Add(new Result(ErrorPlotResultName, new DataTable(ErrorPlotResultName, "Development of errors during gradient descent")));
+      else results[ErrorPlotResultName].Value = new DataTable(ErrorPlotResultName, "Development of errors during gradient descent");
+
+      var plot = results[ErrorPlotResultName].Value as DataTable;
+      if (plot == null) throw new ArgumentException("could not create/access error data table in results collection");
+
+      if (!plot.Rows.ContainsKey("errors")) plot.Rows.Add(new DataRow("errors"));
+      plot.Rows["errors"].Values.Clear();
+      plot.Rows["errors"].VisualProperties.StartIndexZero = true;
+
+      results.Add(new Result(ScatterPlotResultName, "Plot of the projected data", new ScatterPlot(DataResultName, "")));
+      results.Add(new Result(DataResultName, "Projected Data", new DoubleMatrix()));
+    }
+
+    private void Analyze(TSNEStatic<double[]>.TSNEState tsneState) {
+      if (Results == null) return;
+      var results = Results;
+      var plot = results[ErrorPlotResultName].Value as DataTable;
+      if (plot == null) throw new ArgumentException("Could not create/access error data table in results collection.");
+      var errors = plot.Rows["errors"].Values;
+      var c = tsneState.EvaluateError();
+      errors.Add(c);
+      ((IntValue)results[IterationResultName].Value).Value = tsneState.iter;
+      ((DoubleValue)results[ErrorResultName].Value).Value = errors.Last();
+
+      var ndata = Normalize(tsneState.newData);
+      results[DataResultName].Value = new DoubleMatrix(ndata);
+      var splot = results[ScatterPlotResultName].Value as ScatterPlot;
+      FillScatterPlot(ndata, splot);
+    }
+
+    private void FillScatterPlot(double[,] lowDimData, ScatterPlot plot) {
+      foreach (var rowName in dataRowNames.Keys) {
+        if (!plot.Rows.ContainsKey(rowName))
+          plot.Rows.Add(dataRows.ContainsKey(rowName) ? dataRows[rowName] : new ScatterPlotDataRow(rowName, "", new List<Point2D<double>>()));
+        plot.Rows[rowName].Points.Replace(dataRowNames[rowName].Select(i => new Point2D<double>(lowDimData[i, 0], lowDimData[i, 1])));
+      }
+    }
+
+    private static double[,] Normalize(double[,] data) {
+      var max = new double[data.GetLength(1)];
+      var min = new double[data.GetLength(1)];
+      var res = new double[data.GetLength(0), data.GetLength(1)];
+      for (var i = 0; i < max.Length; i++) max[i] = min[i] = data[0, i];
+      for (var i = 0; i < data.GetLength(0); i++)
+        for (var j = 0; j < data.GetLength(1); j++) {
+          var v = data[i, j];
+          max[j] = Math.Max(max[j], v);
+          min[j] = Math.Min(min[j], v);
+        }
+      for (var i = 0; i < data.GetLength(0); i++) {
+        for (var j = 0; j < data.GetLength(1); j++) {
+          res[i, j] = (data[i, j] - (max[j] + min[j]) / 2) / (max[j] - min[j]);
+        }
+      }
+      return res;
+    }
+
+    private static double[][] NormalizeData(IReadOnlyList<double[]> data) {
+      // as in tSNE implementation by van der Maaten
+      var n = data[0].Length;
+      var mean = new double[n];
+      var max = new double[n];
+      var nData = new double[data.Count][];
+      for (var i = 0; i < n; i++) {
+        mean[i] = Enumerable.Range(0, data.Count).Select(x => data[x][i]).Average();
+        max[i] = Enumerable.Range(0, data.Count).Max(x => Math.Abs(data[x][i]));
+      }
+      for (var i = 0; i < data.Count; i++) {
+        nData[i] = new double[n];
+        for (var j = 0; j < n; j++) nData[i][j] = (data[i][j] - mean[j]) / max[j];
+      }
+      return nData;
+    }
+
+    private static Color GetHeatMapColor(int contourNr, int noContours) {
+      var q = (double)contourNr / noContours;  // q in [0,1]
+      var c = q < 0.5 ? Color.FromArgb((int)(q * 2 * 255), 255, 0) : Color.FromArgb(255, (int)((1 - q) * 2 * 255), 0);
+      return c;
+    }
+
+    private static string GetContourName(double value, double min, double max, int noContours) {
+      var size = (max - min) / noContours;
+      var contourNr = (int)((value - min) / size);
+      return GetContourName(contourNr, min, max, noContours);
+    }
+
+    private static string GetContourName(int i, double min, double max, int noContours) {
+      var size = (max - min) / noContours;
+      return "[" + (min + i * size) + ";" + (min + (i + 1) * size) + ")";
+    }
+  }
+}
Index: branches/RBFRegression/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/TSNEStatic.cs
===================================================================
--- branches/RBFRegression/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/TSNEStatic.cs	(revision 14870)
+++ branches/RBFRegression/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/TSNEStatic.cs	(revision 14870)
@@ -0,0 +1,718 @@
+﻿#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.Collections;
+using HeuristicLab.Common;
+using HeuristicLab.Core;
+using HeuristicLab.Optimization;
+using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
+using HeuristicLab.Random;
+
+namespace HeuristicLab.Algorithms.DataAnalysis {
+  [StorableClass]
+  public class TSNEStatic<T> {
+
+    [StorableClass]
+    public sealed class TSNEState : DeepCloneable {
+      // initialized once
+      [Storable]
+      public IDistance<T> distance;
+      [Storable]
+      public IRandom random;
+      [Storable]
+      public double perplexity;
+      [Storable]
+      public bool exact;
+      [Storable]
+      public int noDatapoints;
+      [Storable]
+      public double finalMomentum;
+      [Storable]
+      public int momSwitchIter;
+      [Storable]
+      public int stopLyingIter;
+      [Storable]
+      public double theta;
+      [Storable]
+      public double eta;
+      [Storable]
+      public int newDimensions;
+
+      // for approximate version: sparse representation of similarity/distance matrix
+      [Storable]
+      public double[] valP; // similarity/distance
+      [Storable]
+      public int[] rowP; // row index
+      [Storable]
+      public int[] colP; // col index
+
+      // for exact version: dense representation of distance/similarity matrix
+      [Storable]
+      public double[,] p;
+
+      // mapped data
+      [Storable]
+      public double[,] newData;
+
+      [Storable]
+      public int iter;
+      [Storable]
+      public double currentMomentum;
+
+      // helper variables (updated in each iteration)
+      [Storable]
+      public double[,] gains;
+      [Storable]
+      public double[,] uY;
+      [Storable]
+      public double[,] dY;
+
+      private TSNEState(TSNEState original, Cloner cloner) : base(original, cloner) {
+        this.distance = cloner.Clone(original.distance);
+        this.random = cloner.Clone(original.random);
+        this.perplexity = original.perplexity;
+        this.exact = original.exact;
+        this.noDatapoints = original.noDatapoints;
+        this.finalMomentum = original.finalMomentum;
+        this.momSwitchIter = original.momSwitchIter;
+        this.stopLyingIter = original.stopLyingIter;
+        this.theta = original.theta;
+        this.eta = original.eta;
+        this.newDimensions = original.newDimensions;
+        if(original.valP != null) {
+          this.valP = new double[original.valP.Length];
+          Array.Copy(original.valP, this.valP, this.valP.Length);
+        }
+        if(original.rowP != null) {
+          this.rowP = new int[original.rowP.Length];
+          Array.Copy(original.rowP, this.rowP, this.rowP.Length);
+        }
+        if(original.colP != null) {
+          this.colP = new int[original.colP.Length];
+          Array.Copy(original.colP, this.colP, this.colP.Length);
+        }
+        if(original.p != null) {
+          this.p = new double[original.p.GetLength(0), original.p.GetLength(1)];
+          Array.Copy(original.p, this.p, this.p.Length);
+        }
+        this.newData = new double[original.newData.GetLength(0), original.newData.GetLength(1)];
+        Array.Copy(original.newData, this.newData, this.newData.Length);
+        this.iter = original.iter;
+        this.currentMomentum = original.currentMomentum;
+        this.gains = new double[original.gains.GetLength(0), original.gains.GetLength(1)];
+        Array.Copy(original.gains, this.gains, this.gains.Length);
+        this.uY = new double[original.uY.GetLength(0), original.uY.GetLength(1)];
+        Array.Copy(original.uY, this.uY, this.uY.Length);
+        this.dY = new double[original.dY.GetLength(0), original.dY.GetLength(1)];
+        Array.Copy(original.dY, this.dY, this.dY.Length);
+      }
+
+      public override IDeepCloneable Clone(Cloner cloner) {
+        return new TSNEState(this, cloner);
+      }
+
+      [StorableConstructor]
+      public TSNEState(bool deserializing) { }
+      public TSNEState(T[] data, IDistance<T> distance, IRandom random, int newDimensions, double perplexity, double theta, int stopLyingIter, int momSwitchIter, double momentum, double finalMomentum, double eta) {
+        this.distance = distance;
+        this.random = random;
+        this.newDimensions = newDimensions;
+        this.perplexity = perplexity;
+        this.theta = theta;
+        this.stopLyingIter = stopLyingIter;
+        this.momSwitchIter = momSwitchIter;
+        this.currentMomentum = momentum;
+        this.finalMomentum = finalMomentum;
+        this.eta = eta;
+
+
+        // initialize
+        noDatapoints = data.Length;
+        if(noDatapoints - 1 < 3 * perplexity)
+          throw new ArgumentException("Perplexity too large for the number of data points!");
+
+        exact = Math.Abs(theta) < double.Epsilon;
+        newData = new double[noDatapoints, newDimensions];
+        dY = new double[noDatapoints, newDimensions];
+        uY = new double[noDatapoints, newDimensions];
+        gains = new double[noDatapoints, newDimensions];
+        for(var i = 0; i < noDatapoints; i++)
+          for(var j = 0; j < newDimensions; j++)
+            gains[i, j] = 1.0;
+
+        p = null;
+        rowP = null;
+        colP = null;
+        valP = null;
+
+        //Calculate Similarities
+        if (exact) p = CalculateExactSimilarites(data, distance, perplexity);
+        else CalculateApproximateSimilarities(data, distance, perplexity, out rowP, out colP, out valP);
+
+        // Lie about the P-values (factor is 4 in the MATLAB implementation)
+        if(exact) for(var i = 0; i < noDatapoints; i++) for(var j = 0; j < noDatapoints; j++) p[i, j] *= 12.0;
+        else for(var i = 0; i < rowP[noDatapoints]; i++) valP[i] *= 12.0;
+
+        // Initialize solution (randomly)
+        var rand = new NormalDistributedRandom(random, 0, 1);
+        for(var i = 0; i < noDatapoints; i++)
+          for(var j = 0; j < newDimensions; j++)
+            newData[i, j] = rand.NextDouble() * .0001;
+      }
+
+      public double EvaluateError() {
+        return exact ?
+          EvaluateErrorExact(p, newData, noDatapoints, newDimensions) :
+          EvaluateErrorApproximate(rowP, colP, valP, newData, theta);
+      }
+
+      private static void CalculateApproximateSimilarities(T[] data, IDistance<T> distance, double perplexity, out int[] rowP, out int[] colP, out double[] valP) {
+        // Compute asymmetric pairwise input similarities
+        ComputeGaussianPerplexity(data, distance, out rowP, out colP, out valP, perplexity, (int)(3 * perplexity));
+        // Symmetrize input similarities
+        int[] sRowP, symColP;
+        double[] sValP;
+        SymmetrizeMatrix(rowP, colP, valP, out sRowP, out symColP, out sValP);
+        rowP = sRowP;
+        colP = symColP;
+        valP = sValP;
+        var sumP = .0;
+        for(var i = 0; i < rowP[data.Length]; i++) sumP += valP[i];
+        for(var i = 0; i < rowP[data.Length]; i++) valP[i] /= sumP;
+      }
+
+      private static double[,] CalculateExactSimilarites(T[] data, IDistance<T> distance, double perplexity) {
+        // Compute similarities
+        var p = new double[data.Length, data.Length];
+        ComputeGaussianPerplexity(data, distance, p, perplexity);
+        // Symmetrize input similarities
+        for(var n = 0; n < data.Length; n++) {
+          for(var m = n + 1; m < data.Length; m++) {
+            p[n, m] += p[m, n];
+            p[m, n] = p[n, m];
+          }
+        }
+        var sumP = .0;
+        for(var i = 0; i < data.Length; i++) for(var j = 0; j < data.Length; j++) sumP += p[i, j];
+        for(var i = 0; i < data.Length; i++) for(var j = 0; j < data.Length; j++) p[i, j] /= sumP;
+        return p;
+      }
+
+      private static void ComputeGaussianPerplexity(IReadOnlyList<T> x, IDistance<T> distance, out int[] rowP, out int[] colP, out double[] valP, double perplexity, int k) {
+        if(perplexity > k) throw new ArgumentException("Perplexity should be lower than k!");
+
+        int n = x.Count;
+        // Allocate the memory we need
+        rowP = new int[n + 1];
+        colP = new int[n * k];
+        valP = new double[n * k];
+        var curP = new double[n - 1];
+        rowP[0] = 0;
+        for(var i = 0; i < n; i++) rowP[i + 1] = rowP[i] + k;
+
+        var objX = new List<IndexedItem<T>>();
+        for(var i = 0; i < n; i++) objX.Add(new IndexedItem<T>(i, x[i]));
+
+        // Build ball tree on data set
+        var tree = new VantagePointTree<IndexedItem<T>>(new IndexedItemDistance<T>(distance), objX);
+
+        // Loop over all points to find nearest neighbors
+        for(var i = 0; i < n; i++) {
+          IList<IndexedItem<T>> indices;
+          IList<double> distances;
+
+          // Find nearest neighbors
+          tree.Search(objX[i], k + 1, out indices, out distances);
+
+          // Initialize some variables for binary search
+          var found = false;
+          var beta = 1.0;
+          var minBeta = double.MinValue;
+          var maxBeta = double.MaxValue;
+          const double tol = 1e-5;  
+
+          // Iterate until we found a good perplexity
+          var iter = 0; double sumP = 0;
+          while(!found && iter < 200) {  
+
+            // Compute Gaussian kernel row
+            for(var m = 0; m < k; m++) curP[m] = Math.Exp(-beta * distances[m + 1]); 
+
+            // Compute entropy of current row
+            sumP = double.Epsilon;
+            for(var m = 0; m < k; m++) sumP += curP[m];
+            var h = .0;
+            for(var m = 0; m < k; m++) h += beta * (distances[m + 1] * curP[m]); 
+            h = h / sumP + Math.Log(sumP);
+
+            // Evaluate whether the entropy is within the tolerance level
+            var hdiff = h - Math.Log(perplexity);
+            if(hdiff < tol && -hdiff < tol) {
+              found = true;
+            } else {
+              if(hdiff > 0) {
+                minBeta = beta;
+                if(maxBeta.IsAlmost(double.MaxValue) || maxBeta.IsAlmost(double.MinValue))
+                  beta *= 2.0;
+                else
+                  beta = (beta + maxBeta) / 2.0;
+              } else {
+                maxBeta = beta;
+                if(minBeta.IsAlmost(double.MinValue) || minBeta.IsAlmost(double.MaxValue))
+                  beta /= 2.0;
+                else
+                  beta = (beta + minBeta) / 2.0;
+              }
+            }
+
+            // Update iteration counter
+            iter++;
+          }
+
+          // Row-normalize current row of P and store in matrix
+          for(var m = 0; m < k; m++) curP[m] /= sumP;
+          for(var m = 0; m < k; m++) {
+            colP[rowP[i] + m] = indices[m + 1].Index;
+            valP[rowP[i] + m] = curP[m];
+          }
+        }
+      }
+      private static void ComputeGaussianPerplexity(T[] x, IDistance<T> distance, double[,] p, double perplexity) {
+        // Compute the distance matrix
+        var dd = ComputeDistances(x, distance);
+
+        int n = x.Length;
+        // Compute the Gaussian kernel row by row
+        for(var i = 0; i < n; i++) {
+          // Initialize some variables
+          var found = false;
+          var beta = 1.0;
+          var minBeta = double.MinValue;
+          var maxBeta = double.MaxValue;
+          const double tol = 1e-5;
+          double sumP = 0;
+
+          // Iterate until we found a good perplexity
+          var iter = 0;
+          while(!found && iter < 200) {      // 200 iterations as in tSNE implementation by van der Maarten
+
+            // Compute Gaussian kernel row
+            for(var m = 0; m < n; m++) p[i, m] = Math.Exp(-beta * dd[i][m]);
+            p[i, i] = double.Epsilon;
+
+            // Compute entropy of current row
+            sumP = double.Epsilon;
+            for(var m = 0; m < n; m++) sumP += p[i, m];
+            var h = 0.0;
+            for(var m = 0; m < n; m++) h += beta * (dd[i][m] * p[i, m]);
+            h = h / sumP + Math.Log(sumP);
+
+            // Evaluate whether the entropy is within the tolerance level
+            var hdiff = h - Math.Log(perplexity);
+            if(hdiff < tol && -hdiff < tol) {
+              found = true;
+            } else {
+              if(hdiff > 0) {
+                minBeta = beta;
+                if(maxBeta.IsAlmost(double.MaxValue) || maxBeta.IsAlmost(double.MinValue))
+                  beta *= 2.0;
+                else
+                  beta = (beta + maxBeta) / 2.0;
+              } else {
+                maxBeta = beta;
+                if(minBeta.IsAlmost(double.MinValue) || minBeta.IsAlmost(double.MaxValue))
+                  beta /= 2.0;
+                else
+                  beta = (beta + minBeta) / 2.0;
+              }
+            }
+
+            // Update iteration counter
+            iter++;
+          }
+
+          // Row normalize P
+          for(var m = 0; m < n; m++) p[i, m] /= sumP;
+        }
+      }
+
+      private static double[][] ComputeDistances(T[] x, IDistance<T> distance) {
+        var res = new double[x.Length][];
+        for(int r = 0; r < x.Length; r++) {
+          var rowV = new double[x.Length];
+          // all distances must be symmetric 
+          for(int c = 0; c < r; c++) {
+            rowV[c] = res[c][r];
+          }
+          rowV[r] = 0.0; // distance to self is zero for all distances
+          for(int c = r + 1; c < x.Length; c++) {
+            rowV[c] = distance.Get(x[r], x[c]);
+          }
+          res[r] = rowV;
+        }
+        return res;
+        // return x.Select(m => x.Select(n => distance.Get(m, n)).ToArray()).ToArray();
+      }
+
+      private static double EvaluateErrorExact(double[,] p, double[,] y, int n, int d) {
+        // Compute the squared Euclidean distance matrix
+        var dd = new double[n, n];
+        var q = new double[n, n];
+        ComputeSquaredEuclideanDistance(y, n, d, dd);
+
+        // Compute Q-matrix and normalization sum
+        var sumQ = double.Epsilon;
+        for(var n1 = 0; n1 < n; n1++) {
+          for(var m = 0; m < n; m++) {
+            if(n1 != m) {
+              q[n1, m] = 1 / (1 + dd[n1, m]);
+              sumQ += q[n1, m];
+            } else q[n1, m] = double.Epsilon;
+          }
+        }
+        for(var i = 0; i < n; i++) for(var j = 0; j < n; j++) q[i, j] /= sumQ;
+
+        // Sum t-SNE error
+        var c = .0;
+        for(var i = 0; i < n; i++)
+          for(var j = 0; j < n; j++) {
+            c += p[i, j] * Math.Log((p[i, j] + float.Epsilon) / (q[i, j] + float.Epsilon));
+          }
+        return c;
+      }
+
+      // The mapping of the approximate tSNE looks good but the error curve never changes.
+      private static double EvaluateErrorApproximate(IReadOnlyList<int> rowP, IReadOnlyList<int> colP, IReadOnlyList<double> valP, double[,] y, double theta) {
+        // Get estimate of normalization term
+        var n = y.GetLength(0);
+        var d = y.GetLength(1);
+        var tree = new SpacePartitioningTree(y);
+        var buff = new double[d];
+        double sumQ = 0.0;
+        for(var i = 0; i < n; i++) tree.ComputeNonEdgeForces(i, theta, buff, ref sumQ);
+
+        // Loop over all edges to compute t-SNE error
+        var c = .0;
+        for(var k = 0; k < n; k++) {
+          for(var i = rowP[k]; i < rowP[k + 1]; i++) {
+            var q = .0;
+            for(var j = 0; j < d; j++) buff[j] = y[k, j];
+            for(var j = 0; j < d; j++) buff[j] -= y[colP[i], j];
+            for(var j = 0; j < d; j++) q += buff[j] * buff[j];
+            q = (1.0 / (1.0 + q)) / sumQ;
+            c += valP[i] * Math.Log((valP[i] + float.Epsilon) / (q + float.Epsilon));
+          }
+        }
+        return c;
+      }
+      private static void SymmetrizeMatrix(IReadOnlyList<int> rowP, IReadOnlyList<int> colP, IReadOnlyList<double> valP, out int[] symRowP, out int[] symColP, out double[] symValP) {
+
+        // Count number of elements and row counts of symmetric matrix
+        var n = rowP.Count - 1;
+        var rowCounts = new int[n];
+        for(var j = 0; j < n; j++) {
+          for(var i = rowP[j]; i < rowP[j + 1]; i++) {
+
+            // Check whether element (col_P[i], n) is present
+            var present = false;
+            for(var m = rowP[colP[i]]; m < rowP[colP[i] + 1]; m++) {
+              if(colP[m] == j) present = true;
+            }
+            if(present) rowCounts[j]++;
+            else {
+              rowCounts[j]++;
+              rowCounts[colP[i]]++;
+            }
+          }
+        }
+        var noElem = 0;
+        for(var i = 0; i < n; i++) noElem += rowCounts[i];
+
+        // Allocate memory for symmetrized matrix
+        symRowP = new int[n + 1];
+        symColP = new int[noElem];
+        symValP = new double[noElem];
+
+        // Construct new row indices for symmetric matrix
+        symRowP[0] = 0;
+        for(var i = 0; i < n; i++) symRowP[i + 1] = symRowP[i] + rowCounts[i];
+
+        // Fill the result matrix
+        var offset = new int[n];
+        for(var j = 0; j < n; j++) {
+          for(var i = rowP[j]; i < rowP[j + 1]; i++) {                                  // considering element(n, colP[i])
+
+            // Check whether element (col_P[i], n) is present
+            var present = false;
+            for(var m = rowP[colP[i]]; m < rowP[colP[i] + 1]; m++) {
+              if(colP[m] != j) continue;
+              present = true;
+              if(j > colP[i]) continue; // make sure we do not add elements twice
+              symColP[symRowP[j] + offset[j]] = colP[i];
+              symColP[symRowP[colP[i]] + offset[colP[i]]] = j;
+              symValP[symRowP[j] + offset[j]] = valP[i] + valP[m];
+              symValP[symRowP[colP[i]] + offset[colP[i]]] = valP[i] + valP[m];
+            }
+
+            // If (colP[i], n) is not present, there is no addition involved
+            if(!present) {
+              symColP[symRowP[j] + offset[j]] = colP[i];
+              symColP[symRowP[colP[i]] + offset[colP[i]]] = j;
+              symValP[symRowP[j] + offset[j]] = valP[i];
+              symValP[symRowP[colP[i]] + offset[colP[i]]] = valP[i];
+            }
+
+            // Update offsets
+            if(present && (j > colP[i])) continue;
+            offset[j]++;
+            if(colP[i] != j) offset[colP[i]]++;
+          }
+        }
+
+        for(var i = 0; i < noElem; i++) symValP[i] /= 2.0;
+      }
+    }
+
+    /// <summary>
+    /// Simple interface to tSNE
+    /// </summary>
+    /// <param name="data"></param>
+    /// <param name="distance">The distance function used to differentiate similar from non-similar points, e.g. Euclidean distance.</param>
+    /// <param name="random">Random number generator</param>
+    /// <param name="newDimensions">Dimensionality of projected space (usually 2 for easy visual analysis).</param>
+    /// <param name="perplexity">Perplexity parameter of tSNE. Comparable to k in a k-nearest neighbour algorithm. Recommended value is floor(number of points /3) or lower</param>
+    /// <param name="iterations">Maximum number of iterations for gradient descent.</param>
+    /// <param name="theta">Value describing how much appoximated gradients my differ from exact gradients. Set to 0 for exact calculation and in [0,1] otherwise. CAUTION: exact calculation of forces requires building a non-sparse N*N matrix where N is the number of data points. This may exceed memory limitations.</param>
+    /// <param name="stopLyingIter">Number of iterations after which p is no longer approximated.</param>
+    /// <param name="momSwitchIter">Number of iterations after which the momentum in the gradient descent is switched.</param>
+    /// <param name="momentum">The initial momentum in the gradient descent.</param>
+    /// <param name="finalMomentum">The final momentum in gradient descent (after momentum switch).</param>
+    /// <param name="eta">Gradient descent learning rate.</param>
+    /// <returns></returns>
+    public static double[,] Run(T[] data, IDistance<T> distance, IRandom random,
+      int newDimensions = 2, double perplexity = 25, int iterations = 1000,
+      double theta = 0,
+      int stopLyingIter = 250, int momSwitchIter = 250, double momentum = .5,
+      double finalMomentum = .8, double eta = 200.0
+      ) {
+      var state = CreateState(data, distance, random, newDimensions, perplexity,
+        theta, stopLyingIter, momSwitchIter, momentum, finalMomentum, eta);
+
+      for(int i = 0; i < iterations - 1; i++) {
+        Iterate(state);
+      }
+      return Iterate(state);
+    }
+
+    public static TSNEState CreateState(T[] data, IDistance<T> distance, IRandom random,
+      int newDimensions = 2, double perplexity = 25, double theta = 0,
+      int stopLyingIter = 250, int momSwitchIter = 250, double momentum = .5,
+      double finalMomentum = .8, double eta = 200.0
+      ) {
+      return new TSNEState(data, distance, random, newDimensions, perplexity, theta, stopLyingIter, momSwitchIter, momentum, finalMomentum, eta);
+    }
+
+
+    public static double[,] Iterate(TSNEState state) {
+      if(state.exact)
+        ComputeExactGradient(state.p, state.newData, state.noDatapoints, state.newDimensions, state.dY);
+      else
+        ComputeApproximateGradient(state.rowP, state.colP, state.valP, state.newData, state.noDatapoints, state.newDimensions, state.dY, state.theta);
+
+      // Update gains
+      for(var i = 0; i < state.noDatapoints; i++) {
+        for(var j = 0; j < state.newDimensions; j++) {
+          state.gains[i, j] = Math.Sign(state.dY[i, j]) != Math.Sign(state.uY[i, j])
+            ? state.gains[i, j] + .2  // +0.2 nd *0.8 are used in two separate implementations of tSNE -> seems to be correct
+            : state.gains[i, j] * .8;
+
+          if(state.gains[i, j] < .01) state.gains[i, j] = .01;
+        }
+      }
+
+
+      // Perform gradient update (with momentum and gains)
+      for(var i = 0; i < state.noDatapoints; i++)
+        for(var j = 0; j < state.newDimensions; j++)
+          state.uY[i, j] = state.currentMomentum * state.uY[i, j] - state.eta * state.gains[i, j] * state.dY[i, j];
+
+      for(var i = 0; i < state.noDatapoints; i++)
+        for(var j = 0; j < state.newDimensions; j++)
+          state.newData[i, j] = state.newData[i, j] + state.uY[i, j];
+
+      // Make solution zero-mean
+      ZeroMean(state.newData);
+
+      // Stop lying about the P-values after a while, and switch momentum
+      if(state.iter == state.stopLyingIter) {
+        if(state.exact)
+          for(var i = 0; i < state.noDatapoints; i++)
+            for(var j = 0; j < state.noDatapoints; j++)
+              state.p[i, j] /= 12.0;
+        else
+          for(var i = 0; i < state.rowP[state.noDatapoints]; i++)
+            state.valP[i] /= 12.0;
+      }
+
+      if(state.iter == state.momSwitchIter)
+        state.currentMomentum = state.finalMomentum;
+
+      state.iter++;
+      return state.newData;
+    }
+
+
+
+    private static void ComputeApproximateGradient(int[] rowP, int[] colP, double[] valP, double[,] y, int n, int d, double[,] dC, double theta) {
+      var tree = new SpacePartitioningTree(y);
+      double sumQ = 0.0;
+      var posF = new double[n, d];
+      var negF = new double[n, d];
+      tree.ComputeEdgeForces(rowP, colP, valP, n, posF);
+      var row = new double[d];
+      for(var n1 = 0; n1 < n; n1++) {
+        Array.Clear(row,0,row.Length);
+        tree.ComputeNonEdgeForces(n1, theta, row, ref sumQ);
+        Buffer.BlockCopy(row, 0, negF, (sizeof(double) * n1 * d), d*sizeof(double));
+      }
+
+      // Compute final t-SNE gradient
+      for (var i = 0; i < n; i++)
+        for(var j = 0; j < d; j++) {
+          dC[i, j] = posF[i, j] - negF[i, j] / sumQ;
+        }
+    }
+
+    private static void ComputeExactGradient(double[,] p, double[,] y, int n, int d, double[,] dC) {
+
+      // Make sure the current gradient contains zeros
+      for(var i = 0; i < n; i++) for(var j = 0; j < d; j++) dC[i, j] = 0.0;
+
+      // Compute the squared Euclidean distance matrix
+      var dd = new double[n, n];
+      ComputeSquaredEuclideanDistance(y, n, d, dd);
+
+      // Compute Q-matrix and normalization sum
+      var q = new double[n, n];
+      var sumQ = .0;
+      for(var n1 = 0; n1 < n; n1++) {
+        for(var m = 0; m < n; m++) {
+          if(n1 == m) continue;
+          q[n1, m] = 1 / (1 + dd[n1, m]);
+          sumQ += q[n1, m];
+        }
+      }
+
+      // Perform the computation of the gradient
+      for(var n1 = 0; n1 < n; n1++) {
+        for(var m = 0; m < n; m++) {
+          if(n1 == m) continue;
+          var mult = (p[n1, m] - q[n1, m] / sumQ) * q[n1, m];
+          for(var d1 = 0; d1 < d; d1++) {
+            dC[n1, d1] += (y[n1, d1] - y[m, d1]) * mult;
+          }
+        }
+      }
+    }
+
+    private static void ComputeSquaredEuclideanDistance(double[,] x, int n, int d, double[,] dd) {
+      var dataSums = new double[n];
+      for(var i = 0; i < n; i++) {
+        for(var j = 0; j < d; j++) {
+          dataSums[i] += x[i, j] * x[i, j];
+        }
+      }
+      for(var i = 0; i < n; i++) {
+        for(var m = 0; m < n; m++) {
+          dd[i, m] = dataSums[i] + dataSums[m];
+        }
+      }
+      for(var i = 0; i < n; i++) {
+        dd[i, i] = 0.0;
+        for(var m = i + 1; m < n; m++) {
+          dd[i, m] = 0.0;
+          for(var j = 0; j < d; j++) {
+            dd[i, m] += (x[i, j] - x[m, j]) * (x[i, j] - x[m, j]);
+          }
+          dd[m, i] = dd[i, m];
+        }
+      }
+    }
+
+    private static void ZeroMean(double[,] x) {
+      // Compute data mean
+      var n = x.GetLength(0);
+      var d = x.GetLength(1);
+      var mean = new double[d];
+      for(var i = 0; i < n; i++) {
+        for(var j = 0; j < d; j++) {
+          mean[j] += x[i, j];
+        }
+      }
+      for(var i = 0; i < d; i++) {
+        mean[i] /= n;
+      }
+      // Subtract data mean
+      for(var i = 0; i < n; i++) {
+        for(var j = 0; j < d; j++) {
+          x[i, j] -= mean[j];
+        }
+      }
+    }
+  }
+}
Index: branches/RBFRegression/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/TSNEUtils.cs
===================================================================
--- branches/RBFRegression/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/TSNEUtils.cs	(revision 14870)
+++ branches/RBFRegression/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/TSNEUtils.cs	(revision 14870)
@@ -0,0 +1,80 @@
+﻿#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/>.
+ */
+#endregion
+
+using System;
+using System.Collections.Generic;
+
+namespace HeuristicLab.Algorithms.DataAnalysis {
+  internal static class TSNEUtils {
+    internal static void ForEach<T>(this IEnumerable<T> sequence, Action<int, T> action) {
+      if (sequence == null) return;
+      if (action == null) throw new ArgumentException("Null Action can not be performed");
+      var i = 0;
+      foreach (var item in sequence) {
+        action(i, item);
+        i++;
+      }
+    }
+
+    internal static IList<T> Swap<T>(this IList<T> list, int indexA, int indexB) {
+      var tmp = list[indexA];
+      list[indexA] = list[indexB];
+      list[indexB] = tmp;
+      return list;
+    }
+
+    internal static int Partition<T>(this IList<T> list, int left, int right, int pivotindex, IComparer<T> comparer) {
+      var pivotValue = list[pivotindex];
+      list.Swap(pivotindex, right);
+      var storeIndex = left;
+      for (var i = left; i < right; i++)
+        if (comparer.Compare(list[i], pivotValue) < 0)
+          list.Swap(storeIndex++, i);
+      list.Swap(right, storeIndex);
+      return storeIndex;
+    }
+
+    /// <summary>
+    /// Quick-Sort based partial ascending sorting
+    /// [0,left-1] not changed;
+    /// [left, n] elements are smaller than or equal to list[n]
+    /// [n, right] elements are grater than or equal to list[n]
+    /// [right, list.Count[ not changed;
+    /// </summary>
+    /// <typeparam name="T"></typeparam>
+    /// <param name="list">list that shall be partially sorted</param>
+    /// <param name="left">left index of sorted list part</param>
+    /// <param name="right">right index of sorted list part</param>
+    /// <param name="n">index around which the partial sorting occures</param>
+    /// <param name="comparer">comparer for list elemnts </param>
+    /// <returns></returns>
+    internal static T NthElement<T>(this IList<T> list, int left, int right, int n, IComparer<T> comparer) {
+      while (true) {
+        if (left == right) return list[left];
+        var pivotindex = left + (int)Math.Floor(new System.Random().Next() % (right - (double)left + 1));
+        pivotindex = list.Partition(left, right, pivotindex, comparer);
+        if (n == pivotindex) return list[n];
+        if (n < pivotindex) right = pivotindex - 1;
+        else left = pivotindex + 1;
+      }
+    }
+  }
+}
Index: branches/RBFRegression/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/VantagePointTree.cs
===================================================================
--- branches/RBFRegression/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/VantagePointTree.cs	(revision 14870)
+++ branches/RBFRegression/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/VantagePointTree.cs	(revision 14870)
@@ -0,0 +1,164 @@
+﻿#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.Collections.Generic;
+using System.Linq;
+using HeuristicLab.Collections;
+using HeuristicLab.Random;
+
+namespace HeuristicLab.Algorithms.DataAnalysis {
+  /// <summary>
+  /// Vantage point tree  (or VP tree) is a metric tree that segregates data in a metric space by choosing 
+  /// a position in the space (the "vantage point") and partitioning the data points into two parts: 
+  /// those points that are nearer to the vantage point than a threshold, and those points that are not. 
+  /// By recursively applying this procedure to partition the data into smaller and smaller sets, a tree 
+  /// data structure is created where neighbors in the tree are likely to be neighbors in the space.
+  /// </summary>
+  /// <typeparam name="T"></typeparam>
+  public class VantagePointTree<T> : IVantagePointTree<T> {
+    #region properties
+    private readonly List<T> items;
+    private readonly Node root;
+    private readonly IDistance<T> distance;
+    #endregion
+
+    public VantagePointTree(IDistance<T> distance, IEnumerable<T> items) : base() {
+      root = null;
+      this.distance = distance;
+      this.items = items.Select(x => x).ToList();
+      root = BuildFromPoints(0, this.items.Count);
+    }
+
+    public void Search(T target, int k, out IList<T> results, out IList<double> distances) {
+      var heap = new PriorityQueue<double, IndexedItem<double>>(double.MaxValue, double.MinValue, k);
+      double tau = double.MaxValue;
+      Search(root, target, k, heap, ref tau);
+      var res = new List<T>();
+      var dist = new List<double>();
+      while (heap.Size > 0) {
+        res.Add(items[heap.PeekMinValue().Index]);          // actually max distance
+        dist.Add(heap.PeekMinValue().Value);
+        heap.RemoveMin();
+      }
+      res.Reverse();  
+      dist.Reverse();
+      results = res;
+      distances = dist;
+    }
+
+    private void Search(Node node, T target, int k, PriorityQueue<double, IndexedItem<double>> heap, ref double tau) {
+      if (node == null) return;
+      var dist = distance.Get(items[node.index], target);
+      if (dist < tau) {
+        if (heap.Size == k) heap.RemoveMin();   // remove furthest node from result list (if we already have k results) 
+        heap.Insert(-dist, new IndexedItem<double>(node.index, dist));      
+        if (heap.Size == k) tau = heap.PeekMinValue().Value;
+      }
+      if (node.left == null && node.right == null) return;
+
+      if (dist < node.threshold) {
+        if (dist - tau <= node.threshold) Search(node.left, target, k, heap, ref tau);   // 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, ref tau);  // 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, ref tau);  // 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, ref tau);   // if there can still be neighbors inside the ball, recursively search left child
+      }
+    }
+
+    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() * (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 sealed class Node {
+      public int index;
+      public double threshold;
+      public Node left;
+      public Node right;
+
+      internal Node() {
+        index = 0;
+        threshold = 0;
+        left = null;
+        right = null;
+      }
+    }
+  }
+}
