Index: /branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/Build.cmd
===================================================================
--- /branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/Build.cmd	(revision 14414)
+++ /branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/Build.cmd	(revision 14414)
@@ -0,0 +1,104 @@
+@ECHO OFF
+
+SET CLEANBEFOREBUILD=
+
+SET SELECTED=
+SET CONFIGURATION=
+SET PLATFORM=
+
+IF "%~1"=="" GOTO :prompt_solution
+
+SET SELECTED=%1
+IF NOT EXIST %SELECTED% (
+  ECHO Solution file %SELECTED% could not be found.
+  GOTO :end
+)
+ECHO Building solution %SELECTED% ...
+GOTO :config_selection
+
+:prompt_solution
+SET /A COUNT=0
+FOR /F "tokens=*" %%A IN ('dir /B *.sln') DO (
+  CALL :forloopbody "%%A"
+)
+
+IF "%COUNT%"=="1" (
+  SET SELECTED=%SOLUTIONS.1%
+  ECHO Building %SOLUTIONS.1% as it is the only solution that was found ...
+  GOTO :config_selection
+)
+
+ECHO Found the following solutions:
+FOR /F "tokens=2* delims=.=" %%A IN ('SET SOLUTIONS.') DO ECHO %%A = %%B
+ECHO.
+SET /P SOLUTIONINDEX=Which solution to build? Type the number: 
+
+SET SELECTED=""
+FOR /F "tokens=2* delims=.=" %%A IN ('SET SOLUTIONS.') DO (
+  IF "%%A"=="%SOLUTIONINDEX%" SET SELECTED=%%B
+)
+
+IF %SELECTED%=="" GOTO :eof
+
+:config_selection
+IF "%~2"==""  GOTO :prompt_config
+
+SET CONFIGURATION=%~2
+ECHO Building configuration %CONFIGURATION% ...
+GOTO :platform_selection
+
+:prompt_config
+SET /P CONFIGURATION=Which configuration to build [Release]: 
+IF "%CONFIGURATION%"=="" SET CONFIGURATION=Release
+
+:platform_selection
+IF "%~3"=="" GOTO :prompt_platform
+  
+SET PLATFORM=%~3
+ECHO Building platform %PLATFORM% ...
+GOTO :clean
+
+:prompt_platform
+SET /P PLATFORM=Which platform to build [Any CPU]: 
+IF "%PLATFORM%"=="" SET PLATFORM=Any CPU
+
+:clean
+IF "%~4"=="" GOTO :prompt_clean
+
+SET CLEANBEFOREBUILD=%~4
+GOTO :main
+
+:prompt_clean
+SET /P CLEANBEFOREBUILD=Would you like to clean before building [n]: 
+IF "%CLEANBEFOREBUILD%"=="" SET CLEANBEFOREBUILD=n
+
+:main
+REM First find the path to the msbuild.exe by performing a registry query
+FOR /F "tokens=1,3 delims=	 " %%A IN ('REG QUERY "HKEY_LOCAL_MACHINE\SOFTWARE\Microsoft\MSBuild\ToolsVersions\4.0"') DO (
+  IF "%%A"=="MSBuildToolsPath" SET MSBUILDPATH=%%B)
+
+REM Then execute msbuild to clean and build the solution
+REM Disable that msbuild creates a cache file of the solution
+SET MSBuildUseNoSolutionCache=1
+REM Run msbuild to clean and then build
+IF "%CLEANBEFOREBUILD%" NEQ "n" (
+  ECHO Cleaning ...
+  %MSBUILDPATH%msbuild.exe %SELECTED% /target:Clean /p:Configuration="%CONFIGURATION%",Platform="%PLATFORM%" /m:2 /nologo /verbosity:q /clp:ErrorsOnly
+)
+ECHO Building ...
+%MSBUILDPATH%msbuild.exe %SELECTED% /target:Build /p:Configuration="%CONFIGURATION%",Platform="%PLATFORM%" /m:2 /nologo /verbosity:q /clp:ErrorsOnly
+
+ECHO.
+ECHO DONE.
+
+:end
+
+PAUSE
+
+GOTO :eof
+
+REM This workaround is necessary so that COUNT gets reevaluated
+:forloopbody
+SET /A COUNT+=1
+SET SOLUTIONS.%COUNT%=%1
+GOTO :eof
Index: /branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/HeuristicLab.Algorithms.DataAnalysis-3.4.sln
===================================================================
--- /branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/HeuristicLab.Algorithms.DataAnalysis-3.4.sln	(revision 14414)
+++ /branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/HeuristicLab.Algorithms.DataAnalysis-3.4.sln	(revision 14414)
@@ -0,0 +1,34 @@
+﻿
+Microsoft Visual Studio Solution File, Format Version 12.00
+# Visual Studio 14
+VisualStudioVersion = 14.0.25420.1
+MinimumVisualStudioVersion = 10.0.40219.1
+Project("{FAE04EC0-301F-11D3-BF4B-00C04F79EFBC}") = "HeuristicLab.Algorithms.DataAnalysis-3.4", "HeuristicLab.Algorithms.DataAnalysis-3.4.csproj", "{2E782078-FA81-4B70-B56F-74CE38DAC6C8}"
+EndProject
+Global
+	GlobalSection(SolutionConfigurationPlatforms) = preSolution
+		Debug|Any CPU = Debug|Any CPU
+		Debug|x64 = Debug|x64
+		Debug|x86 = Debug|x86
+		Release|Any CPU = Release|Any CPU
+		Release|x64 = Release|x64
+		Release|x86 = Release|x86
+	EndGlobalSection
+	GlobalSection(ProjectConfigurationPlatforms) = postSolution
+		{2E782078-FA81-4B70-B56F-74CE38DAC6C8}.Debug|Any CPU.ActiveCfg = Debug|Any CPU
+		{2E782078-FA81-4B70-B56F-74CE38DAC6C8}.Debug|Any CPU.Build.0 = Debug|Any CPU
+		{2E782078-FA81-4B70-B56F-74CE38DAC6C8}.Debug|x64.ActiveCfg = Debug|x64
+		{2E782078-FA81-4B70-B56F-74CE38DAC6C8}.Debug|x64.Build.0 = Debug|x64
+		{2E782078-FA81-4B70-B56F-74CE38DAC6C8}.Debug|x86.ActiveCfg = Debug|x86
+		{2E782078-FA81-4B70-B56F-74CE38DAC6C8}.Debug|x86.Build.0 = Debug|x86
+		{2E782078-FA81-4B70-B56F-74CE38DAC6C8}.Release|Any CPU.ActiveCfg = Release|Any CPU
+		{2E782078-FA81-4B70-B56F-74CE38DAC6C8}.Release|Any CPU.Build.0 = Release|Any CPU
+		{2E782078-FA81-4B70-B56F-74CE38DAC6C8}.Release|x64.ActiveCfg = Release|x64
+		{2E782078-FA81-4B70-B56F-74CE38DAC6C8}.Release|x64.Build.0 = Release|x64
+		{2E782078-FA81-4B70-B56F-74CE38DAC6C8}.Release|x86.ActiveCfg = Release|x86
+		{2E782078-FA81-4B70-B56F-74CE38DAC6C8}.Release|x86.Build.0 = Release|x86
+	EndGlobalSection
+	GlobalSection(SolutionProperties) = preSolution
+		HideSolutionNode = FALSE
+	EndGlobalSection
+EndGlobal
Index: /branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/PreBuildEvent.cmd
===================================================================
--- /branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/PreBuildEvent.cmd	(revision 14414)
+++ /branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/PreBuildEvent.cmd	(revision 14414)
@@ -0,0 +1,11 @@
+IF EXIST "%ProjectDir%\Properties\AssemblyInfo.cs.frame" SubWCRev "%ProjectDir%\" "%ProjectDir%\Properties\AssemblyInfo.cs.frame" "%ProjectDir%\Properties\AssemblyInfo.cs"
+IF %ERRORLEVEL% NEQ 0 GOTO Error_Handling
+IF EXIST "%ProjectDir%\Plugin.cs.frame" SubWCRev "%ProjectDir%\" "%ProjectDir%\Plugin.cs.frame" "%ProjectDir%\Plugin.cs"
+IF %ERRORLEVEL% NEQ 0 GOTO Error_Handling
+GOTO Done
+
+:Error_Handling
+ECHO There was an error while running subwcrev. Please verify that the *.cs.frame files have been correctly converted to *.cs files, otherwise HeuristicLab won't build. 
+exit 0
+
+:Done
Index: /branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/Cell.cs
===================================================================
--- /branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/Cell.cs	(revision 14414)
+++ /branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/Cell.cs	(revision 14414)
@@ -0,0 +1,117 @@
+﻿#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.Common;
+using HeuristicLab.Core;
+using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
+
+namespace HeuristicLab.Algorithms.DataAnalysis {
+  [StorableClass]
+  public class Cell : Item, ICell {
+    #region properties
+    [Storable]
+    private uint dimension;
+    [Storable]
+    private double[] corner;
+    [Storable]
+    private double[] width;
+    #endregion
+
+    #region HLConstructors & Cloning
+    [StorableConstructor]
+    protected Cell(bool deserializing) : base(deserializing) { }
+    protected Cell(Cell original, Cloner cloner) : base(original, cloner) {
+      dimension = original.dimension;
+      corner = original.corner.Select(x => x).ToArray();
+      width = original.width.Select(x => x).ToArray();
+    }
+    public override IDeepCloneable Clone(Cloner cloner) { return new Cell(this, cloner); }
+
+    public Cell(uint inpDimension) {
+      dimension = inpDimension;
+      corner = new double[dimension];
+      width = new double[dimension];
+    }
+    public Cell(uint inpDimension, IEnumerable<double> inpCorner, IEnumerable<double> inpWidth) {
+      dimension = inpDimension;
+      corner = inpCorner.Select(x => x).ToArray();
+      width = inpWidth.Select(x => x).ToArray();
+    }
+    #endregion
+
+    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/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/DataPoint.cs
===================================================================
--- /branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/DataPoint.cs	(revision 14414)
+++ /branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/DataPoint.cs	(revision 14414)
@@ -0,0 +1,94 @@
+﻿#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 HeuristicLab.Common;
+using HeuristicLab.Core;
+using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
+
+namespace HeuristicLab.Algorithms.DataAnalysis {
+  [StorableClass]
+  public class DataPoint<T> : Item, IDataPoint<T> where T : class, IDeepCloneable {
+    #region properties
+    [Storable]
+    private int ind;
+    [Storable]
+    public T X { get; private set; }
+
+    public int Index
+    {
+      get { return ind; }
+    }
+    #endregion
+
+    #region HLConstructors & Cloning
+    [StorableConstructor]
+    protected DataPoint(bool deserializing) : base(deserializing) { }
+    protected DataPoint(DataPoint<T> original, Cloner cloner)
+      : base(original, cloner) {
+      ind = original.ind;
+      X = cloner.Clone(original.X);
+    }
+    public override IDeepCloneable Clone(Cloner cloner) { return new DataPoint<T>(this, cloner); }
+    public DataPoint() : this(1, null) { }
+    public DataPoint(int ind, T x) {
+      this.ind = ind;
+      X = (T)x.Clone();
+    }
+    #endregion
+
+
+
+  }
+}
Index: /branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/Distances/DataPointDistance.cs
===================================================================
--- /branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/Distances/DataPointDistance.cs	(revision 14414)
+++ /branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/Distances/DataPointDistance.cs	(revision 14414)
@@ -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 HeuristicLab.Algorithms.DataAnalysis.Distances;
+using HeuristicLab.Common;
+using HeuristicLab.Core;
+using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
+
+namespace HeuristicLab.Algorithms.DataAnalysis {
+  [StorableClass]
+  [Item("DataPointDistance", "A distance wrapper for DataPoints")]
+  internal class DataPointDistance<T> : DistanceBase<IDataPoint<T>> where T : class {
+
+    #region properties
+    [Storable]
+    private IDistance<T> dist;
+    #endregion
+
+    #region HLConstructors & Cloning
+    [StorableConstructor]
+    protected DataPointDistance(bool deserializing) : base(deserializing) { }
+    protected DataPointDistance(DataPointDistance<T> original, Cloner cloner) : base(original, cloner) {
+      dist = cloner.Clone(original.dist);
+    }
+    public override IDeepCloneable Clone(Cloner cloner) { return new DataPointDistance<T>(this, cloner); }
+    public DataPointDistance(IDistance<T> distance) {
+      dist = distance;
+    }
+    #endregion
+
+    public override double Get(IDataPoint<T> a, IDataPoint<T> b) {
+      return dist.Get(a.X, b.X);
+    }
+
+    public override double Get(IDataPoint<T> a, IDataPoint<T> b, double threshold) {
+      return dist.Get(a.X, b.X, threshold);
+    }
+  }
+}
+
Index: /branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/Distances/DistanceBase.cs
===================================================================
--- /branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/Distances/DistanceBase.cs	(revision 14414)
+++ /branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/Distances/DistanceBase.cs	(revision 14414)
@@ -0,0 +1,59 @@
+﻿#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.Distances {
+  [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 abstract double Get(T a, T b, double threshold);
+
+    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/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/Distances/EuclidianDistance.cs
===================================================================
--- /branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/Distances/EuclidianDistance.cs	(revision 14414)
+++ /branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/Distances/EuclidianDistance.cs	(revision 14414)
@@ -0,0 +1,69 @@
+﻿#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.Algorithms.DataAnalysis.Distances;
+using HeuristicLab.Common;
+using HeuristicLab.Core;
+using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
+
+namespace HeuristicLab.Algorithms.DataAnalysis {
+  [StorableClass]
+  [Item("EuclidianNorm", "A norm function that uses Euclidean distance")]
+  public class EuclidianDistance : DistanceBase<IEnumerable<double>> {
+
+    #region HLConstructors & Boilerplate
+    [StorableConstructor]
+    protected EuclidianDistance(bool deserializing) : base(deserializing) { }
+    protected EuclidianDistance(EuclidianDistance original, Cloner cloner) : base(original, cloner) { }
+    public override IDeepCloneable Clone(Cloner cloner) { return new EuclidianDistance(this, cloner); }
+    public EuclidianDistance() { }
+    #endregion
+
+    #region statics
+    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 static double GetDistance(IEnumerable<double> a, IEnumerable<double> b, double threshold) {
+      double sum = 0;
+      var it1 = a.GetEnumerator();
+      var it2 = b.GetEnumerator();
+      while (sum > threshold * threshold && it1.MoveNext() && it2.MoveNext()) {
+        var d = it1.Current - it2.Current;
+        sum += d * d;
+      }
+      it1.Dispose();
+      it2.Dispose();
+      return sum;
+    }
+    #endregion
+
+    public override double Get(IEnumerable<double> a, IEnumerable<double> b) {
+      return GetDistance(a.ToArray(), b.ToArray());
+    }
+    public override double Get(IEnumerable<double> a, IEnumerable<double> b, double threshold) {
+      return GetDistance(a, b, threshold);
+    }
+  }
+}
Index: /branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/PriorityQueue.cs
===================================================================
--- /branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/PriorityQueue.cs	(revision 14414)
+++ /branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/PriorityQueue.cs	(revision 14414)
@@ -0,0 +1,203 @@
+﻿#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 {
+
+  // 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> : IHeap<TK, TV> where TK : IComparable {
+    private class KNElement {
+      public TK Key { get; set; }
+      public TV Value { get; set; }
+    }
+
+    private 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 KNElement[] data;    // alligned 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/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/SPtree.cs
===================================================================
--- /branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/SPtree.cs	(revision 14414)
+++ /branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/SPtree.cs	(revision 14414)
@@ -0,0 +1,329 @@
+﻿#region License Information
+/* HeuristicLab
+ * Copyright (C) 2002-2016 Heuristic and Evolutionary Algorithms Laboratory (HEAL)
+ *
+ * This file is part of HeuristicLab.
+ *
+ * HeuristicLab is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * HeuristicLab is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with HeuristicLab. If not, see <http://www.gnu.org/licenses/>.
+ */
+
+//Code is based on an implementation from Laurens van der Maaten
+
+/*
+*
+* Copyright (c) 2014, Laurens van der Maaten (Delft University of Technology)
+* All rights reserved.
+*
+* Redistribution and use in source and binary forms, with or without
+* modification, are permitted provided that the following conditions are met:
+* 1. Redistributions of source code must retain the above copyright
+*    notice, this list of conditions and the following disclaimer.
+* 2. Redistributions in binary form must reproduce the above copyright
+*    notice, this list of conditions and the following disclaimer in the
+*    documentation and/or other materials provided with the distribution.
+* 3. All advertising materials mentioning features or use of this software
+*    must display the following acknowledgement:
+*    This product includes software developed by the Delft University of Technology.
+* 4. Neither the name of the Delft University of Technology nor the names of
+*    its contributors may be used to endorse or promote products derived from
+*    this software without specific prior written permission.
+*
+* THIS SOFTWARE IS PROVIDED BY LAURENS VAN DER MAATEN ''AS IS'' AND ANY EXPRESS
+* OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
+* OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO
+* EVENT SHALL LAURENS VAN DER MAATEN BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+* SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+* PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR
+* BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+* CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING
+* IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY
+* OF SUCH DAMAGE.
+*
+*/
+#endregion
+
+using System;
+using System.Collections.Generic;
+using System.Linq;
+using HeuristicLab.Common;
+using HeuristicLab.Core;
+using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
+
+namespace HeuristicLab.Algorithms.DataAnalysis {
+  [StorableClass]
+  public class SPTree : Item, ISPTree {
+    private const uint QT_NODE_CAPACITY = 1;
+
+    #region properties
+    [Storable]
+    private double[] buff;
+    [Storable]
+    private ISPTree parent;
+    [Storable]
+    private int dimension;
+    [Storable]
+    private bool isLeaf;
+    [Storable]
+    private uint size;
+    [Storable]
+    private uint cumulativeSize;
+
+    // Axis-aligned bounding box stored as a center with half-dimensions to represent the boundaries of this quad tree
+    [Storable]
+    private ICell boundary;
+
+    // Indices in this space-partitioning tree node, corresponding center-of-mass, and list of all children
+    [Storable]
+    public double[,] Data { private get; set; }
+    [Storable]
+    private double[] centerOfMass;
+    [Storable]
+    private int[] index = new int[QT_NODE_CAPACITY];
+
+    // Children
+    private ISPTree[] children;
+    [Storable]
+    private uint noChildren;
+    #endregion
+
+    #region HLConstructors & Cloning
+    [StorableConstructor]
+    protected SPTree(bool deserializing) : base(deserializing) { }
+    protected SPTree(SPTree original, Cloner cloner)
+      : base(original, cloner) {
+      buff = original.buff.Select(x => x).ToArray();
+      parent = cloner.Clone(original.parent);
+      dimension = original.dimension;
+      isLeaf = original.isLeaf;
+      size = original.size;
+      cumulativeSize = original.cumulativeSize;
+      boundary = cloner.Clone(original.boundary);
+      if (original.Data != null) {
+        Data = new double[size, dimension];
+        Buffer.BlockCopy(original.Data, 0, Data, 0, original.Data.Length * sizeof(double));
+      }
+      centerOfMass = original.centerOfMass.Select(x => x).ToArray();
+      index = original.index.Select(x => x).ToArray();
+      children = original.children.Select(cloner.Clone).ToArray();
+      noChildren = original.noChildren;
+    }
+    public override IDeepCloneable Clone(Cloner cloner) { return new SPTree(this, cloner); }
+
+    public SPTree(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] - maxY[i]) + 1e-5;
+      Init(null, inpData, meanY, width);
+      Fill(n);
+    }
+
+    public SPTree(double[,] inpData, IEnumerable<double> impCorner, IEnumerable<double> impWith) {
+      Init(null, inpData, impCorner, impWith);
+    }
+    public SPTree(ISPTree parent, double[,] inpData, IEnumerable<double> impCorner, IEnumerable<double> impWith) {
+      Init(parent, inpData, impCorner, impWith);
+    }
+    #endregion
+
+    public ISPTree 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, (int)(sizeof(double) * dimension * newIndex), point, 0, (int)(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 SPTree(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] = int.MaxValue;
+      }
+
+      // 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, (int)(sizeof(double) * dimension * n), row, 0, (int)(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, 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[0] += 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, sumQ);
+      }
+    }
+
+    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
+          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(ISPTree 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 ISPTree[noChildren];
+      centerOfMass = new double[dimension];
+      buff = new double[dimension];
+
+    }
+    #endregion
+
+  }
+}
Index: /branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/TSNE.cs
===================================================================
--- /branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/TSNE.cs	(revision 14414)
+++ /branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/TSNE.cs	(revision 14414)
@@ -0,0 +1,588 @@
+﻿#region License Information
+/* HeuristicLab
+ * Copyright (C) 2002-2016 Heuristic and Evolutionary Algorithms Laboratory (HEAL)
+ *
+ * This file is part of HeuristicLab.
+ *
+ * HeuristicLab is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * HeuristicLab is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with HeuristicLab. If not, see <http://www.gnu.org/licenses/>.
+ */
+
+//Code is based on an implementation from Laurens van der Maaten
+
+/*
+*
+* Copyright (c) 2014, Laurens van der Maaten (Delft University of Technology)
+* All rights reserved.
+*
+* Redistribution and use in source and binary forms, with or without
+* modification, are permitted provided that the following conditions are met:
+* 1. Redistributions of source code must retain the above copyright
+*    notice, this list of conditions and the following disclaimer.
+* 2. Redistributions in binary form must reproduce the above copyright
+*    notice, this list of conditions and the following disclaimer in the
+*    documentation and/or other materials provided with the distribution.
+* 3. All advertising materials mentioning features or use of this software
+*    must display the following acknowledgement:
+*    This product includes software developed by the Delft University of Technology.
+* 4. Neither the name of the Delft University of Technology nor the names of
+*    its contributors may be used to endorse or promote products derived from
+*    this software without specific prior written permission.
+*
+* THIS SOFTWARE IS PROVIDED BY LAURENS VAN DER MAATEN ''AS IS'' AND ANY EXPRESS
+* OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
+* OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO
+* EVENT SHALL LAURENS VAN DER MAATEN BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+* SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+* PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR
+* BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+* CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING
+* IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY
+* OF SUCH DAMAGE.
+*
+*/
+#endregion
+
+using System;
+using System.Collections.Generic;
+using System.Linq;
+using HeuristicLab.Common;
+using HeuristicLab.Core;
+using HeuristicLab.Data;
+using HeuristicLab.Optimization;
+using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
+using HeuristicLab.Random;
+
+namespace HeuristicLab.Algorithms.DataAnalysis {
+  [StorableClass]
+  public class TSNE<T> : Item, ITSNE<T> where T : class, IDeepCloneable {
+
+    private const string IterationResultName = "Iteration";
+    private const string ErrorResultName = "Error";
+
+    #region Properties
+    [Storable]
+    private IDistance<T> distance;
+    [Storable]
+    private int maxIter;
+    [Storable]
+    private int stopLyingIter;
+    [Storable]
+    private int momSwitchIter;
+    [Storable]
+    double momentum;
+    [Storable]
+    private double finalMomentum;
+    [Storable]
+    private double eta;
+    [Storable]
+    private IRandom random;
+    [Storable]
+    private ResultCollection results;
+    #endregion
+
+    #region HLConstructors & Cloning
+    [StorableConstructor]
+    protected TSNE(bool deserializing) : base(deserializing) { }
+    protected TSNE(TSNE<T> original, Cloner cloner) : base(original, cloner) {
+      distance = cloner.Clone(original.distance);
+      maxIter = original.maxIter;
+      stopLyingIter = original.stopLyingIter;
+      momSwitchIter = original.momSwitchIter;
+      momentum = original.momentum;
+      finalMomentum = original.finalMomentum;
+      eta = original.eta;
+      random = cloner.Clone(random);
+      results = cloner.Clone(results);
+    }
+    public override IDeepCloneable Clone(Cloner cloner) { return new TSNE<T>(this, cloner); }
+    public TSNE(IDistance<T> distance, IRandom random, ResultCollection results = null, int maxIter = 1000, int stopLyingIter = 250, int momSwitchIter = 250, double momentum = .5, double finalMomentum = .8, double eta = 200.0) {
+      this.distance = distance;
+      this.maxIter = maxIter;
+      this.stopLyingIter = stopLyingIter;
+      this.momSwitchIter = momSwitchIter;
+      this.momentum = momentum;
+      this.finalMomentum = finalMomentum;
+      this.eta = eta;
+      this.random = random;
+      this.results = results;
+    }
+    #endregion
+
+    public double[,] Run(T[] data, int newDimensions, double perplexity, double theta) {
+      var currentMomentum = momentum;
+      var noDatapoints = data.Length;
+      if (noDatapoints - 1 < 3 * perplexity) throw new ArgumentException("Perplexity too large for the number of data points!");
+
+      if (results != null) {
+        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;
+      }
+
+      // Determine whether we are using an exact algorithm
+      var exact = Math.Abs(theta) < double.Epsilon;
+      var newData = new double[noDatapoints, newDimensions];
+      var dY = new double[noDatapoints, newDimensions];
+      var uY = new double[noDatapoints, newDimensions];
+      var gains = new double[noDatapoints, newDimensions];
+      for (var i = 0; i < noDatapoints; i++) for (var j = 0; j < newDimensions; j++) gains[i, j] = 1.0;
+
+      // Compute input similarities for exact t-SNE
+      double[,] p = null;
+      int[] rowP = null;
+      int[] colP = null;
+      double[] valP = null;
+      if (exact) {
+        // Compute similarities
+        p = new double[noDatapoints, noDatapoints];
+        ComputeGaussianPerplexity(data, noDatapoints, p, perplexity);
+        // Symmetrize input similarities
+        for (var n = 0; n < noDatapoints; n++) {
+          for (var m = n + 1; m < noDatapoints; m++) {
+            p[n, m] += p[m, n];
+            p[m, n] = p[n, m];
+          }
+        }
+        var sumP = .0;
+        for (var i = 0; i < noDatapoints; i++) for (var j = 0; j < noDatapoints; j++) sumP += p[i, j];
+        for (var i = 0; i < noDatapoints; i++) for (var j = 0; j < noDatapoints; j++) p[i, j] /= sumP;
+      } // Compute input similarities for approximate t-SNE
+      else {
+        // Compute asymmetric pairwise input similarities
+        ComputeGaussianPerplexity(data, noDatapoints, 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[noDatapoints]; i++) sumP += valP[i];
+        for (var i = 0; i < rowP[noDatapoints]; i++) valP[i] /= sumP;
+      }
+
+      // Lie about the P-values
+      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;
+      }
+
+      var rand = new NormalDistributedRandom(random, 0, 1);
+      // Initialize solution (randomly)
+      for (var i = 0; i < noDatapoints; i++)
+        for (var j = 0; j < newDimensions; j++)
+          newData[i, j] = rand.NextDouble() * .0001;
+
+      // Perform main training loop
+      for (var iter = 0; iter < maxIter; iter++) {
+
+        // Compute (approximate) gradient
+        if (exact) ComputeExactGradient(p, newData, noDatapoints, newDimensions, dY);
+        else ComputeGradient(rowP, colP, valP, newData, noDatapoints, newDimensions, dY, theta);
+
+        // Update gains
+        for (var i = 0; i < noDatapoints; i++) for (var j = 0; j < newDimensions; j++)
+            gains[i, j] = Math.Sign(dY[i, j]) != Math.Sign(uY[i, j]) ? gains[i, j] + .2 : gains[i, j] * .8;
+        for (var i = 0; i < noDatapoints; i++) for (var j = 0; j < newDimensions; j++) if (gains[i, j] < .01) gains[i, j] = .01;
+
+        // Perform gradient update (with momentum and gains)
+        for (var i = 0; i < noDatapoints; i++) for (var j = 0; j < newDimensions; j++) uY[i, j] = currentMomentum * uY[i, j] - eta * gains[i, j] * dY[i, j];
+        for (var i = 0; i < noDatapoints; i++) for (var j = 0; j < newDimensions; j++) newData[i, j] = newData[i, j] + uY[i, j];
+
+        // Make solution zero-mean
+        ZeroMean(newData);
+
+        // Stop lying about the P-values after a while, and switch momentum
+        if (iter == stopLyingIter) {
+          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;
+          }
+        }
+        if (iter == momSwitchIter) currentMomentum = finalMomentum;
+
+        if (results == null) continue;
+        var errors = new List<double>();
+        // Print out progress
+        var c = exact
+          ? EvaluateError(p, newData, noDatapoints, newDimensions)
+          : EvaluateError(rowP, colP, valP, newData, theta);
+        errors.Add(c);
+        ((IntValue)results[IterationResultName].Value).Value = iter + 1;
+        ((DoubleValue)results[ErrorResultName].Value).Value = errors.Last();
+      }
+      return newData;
+    }
+    public static double[,] Run<TR>(TR[] data, int newDimensions, double perplexity, double theta, IDistance<TR> distance, IRandom random) where TR : class, IDeepCloneable {
+      return new TSNE<TR>(distance, random).Run(data, newDimensions, perplexity, theta);
+    }
+
+    #region helpers
+    private void ComputeGaussianPerplexity(IReadOnlyList<T> x, int n, 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!");
+
+      // 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;
+
+      // Build ball tree on data set
+      var tree = new VPTree<IDataPoint<T>>(new DataPointDistance<T>(distance));
+      var objX = new List<IDataPoint<T>>();
+      for (var i = 0; i < n; i++) objX.Add(new DataPoint<T>(i, x[i]));
+      tree.Create(objX);
+
+      // Loop over all points to find nearest neighbors
+      var indices = new List<IDataPoint<T>>();
+      var distances = new List<double>();
+      for (var i = 0; i < n; i++) {
+
+        // Find nearest neighbors
+        indices.Clear();
+        distances.Clear();
+        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.MaxValue;
+        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 void ComputeGaussianPerplexity(T[] x, int n, double[,] p, double perplexity) {
+      // Compute the squared Euclidean distance matrix
+      var dd = ComputeDistances(x);
+      // 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.MaxValue;
+        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) {
+
+          // 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 double[][] ComputeDistances(T[] x) {
+      return x.Select(m => x.Select(n => distance.Get(m, n)).ToArray()).ToArray();
+    }
+    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 ComputeGradient(int[] rowP, int[] colP, double[] valP, double[,] y, int n, int d, double[,] dC, double theta) {
+      var tree = new SPTree(y);
+      double[] sumQ = { 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 (int n1 = 0; n1 < n; n1++) {
+        Buffer.BlockCopy(negF, (sizeof(double) * n1 * d), row, 0, d);
+        tree.ComputeNonEdgeForces(n1, theta, row, sumQ);
+      }
+
+      // 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[0];
+        }
+    }
+    private static double EvaluateError(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;
+    }
+    private static double EvaluateError(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 SPTree(y);
+      var buff = new double[d];
+      double[] sumQ = { 0 };
+      for (var i = 0; i < n; i++) tree.ComputeNonEdgeForces(i, theta, buff, 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[0];
+          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]]++;
+        }
+      }
+
+      // Divide the result by two
+      for (var i = 0; i < noElem; i++) symValP[i] /= 2.0;
+    }
+    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];
+        }
+      }
+    }
+    #endregion
+  }
+}
Index: /branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/TSNEAnalysis.cs
===================================================================
--- /branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/TSNEAnalysis.cs	(revision 14414)
+++ /branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/TSNEAnalysis.cs	(revision 14414)
@@ -0,0 +1,230 @@
+#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.Linq;
+using HeuristicLab.Analysis;
+using HeuristicLab.Common;
+using HeuristicLab.Core;
+using HeuristicLab.Data;
+using HeuristicLab.Encodings.RealVectorEncoding;
+using HeuristicLab.Optimization;
+using HeuristicLab.Parameters;
+using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
+using HeuristicLab.Problems.DataAnalysis;
+using HeuristicLab.Random;
+
+namespace HeuristicLab.Algorithms.DataAnalysis {
+  /// <summary>
+  /// Linear regression data analysis algorithm.
+  /// </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 TSNEAnalysis : FixedDataAnalysisAlgorithm<IRegressionProblem> {
+
+    #region Resultnames
+    private const string ScatterPlotResultName = "Scatterplot";
+    private const string DataResultName = "Projected Data";
+    #endregion
+
+    #region Parameternames
+    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";
+    #endregion
+
+    #region Parameterproperties
+    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<RealVector>> DistanceParameter
+    {
+      get { return Parameters[DistanceParameterName] as IValueParameter<IDistance<RealVector>>; }
+    }
+    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>; }
+    }
+    #endregion
+
+    #region  Properties
+    public IDistance<RealVector> Distance
+    {
+      get { return DistanceParameter.Value; }
+    }
+    public double Perplexity
+    {
+      get { return PerplexityParameter.Value.Value; }
+    }
+    public double Theta
+    {
+      get { return ThetaParameter.Value.Value; }
+    }
+    public int NewDimensions
+    {
+      get { return NewDimensionsParameter.Value.Value; }
+    }
+    public int MaxIterations
+    {
+      get { return MaxIterationsParameter.Value.Value; }
+    }
+    public int StopLyingIteration
+    {
+      get { return StopLyingIterationParameter.Value.Value; }
+    }
+    public int MomentumSwitchIteration
+    {
+      get { return MomentumSwitchIterationParameter.Value.Value; }
+    }
+    public double InitialMomentum
+    {
+      get { return InitialMomentumParameter.Value.Value; }
+    }
+    public double FinalMomentum
+    {
+      get { return FinalMomentumParameter.Value.Value; }
+    }
+    public double Eta
+    {
+      get { return EtaParameter.Value.Value; }
+    }
+    public bool SetSeedRandomly
+    {
+      get { return SetSeedRandomlyParameter.Value.Value; }
+    }
+    public uint Seed
+    {
+      get { return (uint)SeedParameter.Value.Value; }
+    }
+    #endregion
+
+    #region Constructors & Cloning
+    [StorableConstructor]
+    private TSNEAnalysis(bool deserializing) : base(deserializing) { }
+    private TSNEAnalysis(TSNEAnalysis original, Cloner cloner) : base(original, cloner) { }
+    public override IDeepCloneable Clone(Cloner cloner) { return new TSNEAnalysis(this, cloner); }
+    public TSNEAnalysis() {
+      Problem = new RegressionProblem();
+      Parameters.Add(new ValueParameter<IDistance<RealVector>>(DistanceParameterName, "The distance function used to differentiate similar from non-similar points", new EuclidianDistance()));
+      Parameters.Add(new FixedValueParameter<DoubleValue>(PerplexityParameterName, "Perplexity-Parameter of TSNE. Comparable to k in a k-nearest neighbour algorithm", 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", new DoubleValue(0.1)));
+      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(250)));
+      Parameters.Add(new FixedValueParameter<IntValue>(MomentumSwitchIterationParameterName, "Number of iterations after which the momentum in the gradient descent is switched", new IntValue(250)));
+      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(200)));
+      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)));
+    }
+    #endregion
+
+    protected override void Run() {
+      var lowDimData = new DoubleMatrix(GetProjectedData(Problem.ProblemData));
+      Results.Add(new Result(ScatterPlotResultName, "Plot of the projected data", CreateScatterPlot(lowDimData, Problem.ProblemData)));
+      Results.Add(new Result(DataResultName, "Projected Data", lowDimData));
+    }
+
+    private ScatterPlot CreateScatterPlot(DoubleMatrix lowDimData, IDataAnalysisProblemData problemData) {
+      var plot = new ScatterPlot(DataResultName, "");
+      Normalize(lowDimData);
+      plot.Rows.Add(new ScatterPlotDataRow("Training", "Points of the training set", problemData.TrainingIndices.Select(i => new Point2D<double>(lowDimData[i, 0], lowDimData[i, 1]))));
+      plot.Rows.Add(new ScatterPlotDataRow("Test", "Points of the test set", problemData.TestIndices.Select(i => new Point2D<double>(lowDimData[i, 0], lowDimData[i, 1]))));
+      return plot;
+    }
+
+    private double[,] GetProjectedData(IDataAnalysisProblemData problemData) {
+      var random = SetSeedRandomly ? new MersenneTwister() : new MersenneTwister(Seed);
+      var tsne = new TSNE<RealVector>(Distance, random, Results, MaxIterations, StopLyingIteration, MomentumSwitchIteration, InitialMomentum, FinalMomentum, Eta);
+      var dataset = problemData.Dataset;
+      var allowedInputVariables = problemData.AllowedInputVariables.ToArray();
+      var data = new RealVector[dataset.Rows];
+      for (var row = 0; row < dataset.Rows; row++) data[row] = new RealVector(allowedInputVariables.Select(col => dataset.GetDoubleValue(col, row)).ToArray());
+      return tsne.Run(data, NewDimensions, Perplexity, Theta);
+    }
+
+    private static void Normalize(DoubleMatrix data) {
+      var max = new double[data.Columns];
+      var min = new double[data.Columns];
+      for (var i = 0; i < max.Length; i++) max[i] = min[i] = data[0, i];
+      for (var i = 0; i < data.Rows; i++)
+        for (var j = 0; j < data.Columns; 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.Rows; i++) {
+        for (var j = 0; j < data.Columns; j++) {
+          data[i, j] = (data[i, j] - (max[j] + min[j]) / 2) / (max[j] - min[j]);
+        }
+      }
+
+    }
+  }
+}
Index: /branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/TSNEInterfaces/ICell.cs
===================================================================
--- /branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/TSNEInterfaces/ICell.cs	(revision 14414)
+++ /branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/TSNEInterfaces/ICell.cs	(revision 14414)
@@ -0,0 +1,32 @@
+﻿#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 ICell : IDeepCloneable {
+    double GetCorner(int d);
+    double GetWidth(int d);
+    void SetCorner(int d, double val);
+    void SetWidth(int d, double val);
+    bool ContainsPoint(double[] point);
+  }
+}
Index: /branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/TSNEInterfaces/IDataPoint.cs
===================================================================
--- /branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/TSNEInterfaces/IDataPoint.cs	(revision 14414)
+++ /branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/TSNEInterfaces/IDataPoint.cs	(revision 14414)
@@ -0,0 +1,30 @@
+﻿#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 IDataPoint<out T> : IDeepCloneable {
+    int Index { get; }
+    T X { get; }
+
+  }
+}
Index: /branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/TSNEInterfaces/IDistance.cs
===================================================================
--- /branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/TSNEInterfaces/IDistance.cs	(revision 14414)
+++ /branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/TSNEInterfaces/IDistance.cs	(revision 14414)
@@ -0,0 +1,55 @@
+﻿#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 aribitrary Vectors. The distance value d is 
+    /// 1.) positive  d(x,y)>=0
+    /// 2.) symmetric d(x,y) = d(y,x)
+    /// 3.) zero-reflexive d(x,x) =0;
+    /// </summary>
+    /// <param name="a">an array representing point x</param>
+    /// <param name="b">>an array representing point y</param>
+    /// <returns>d(x,y)</returns>
+    double Get(T a, T b);
+
+    /// <summary>
+    /// Calculates the correct kernel measure if it is smaller than threshold, but any value greater then threshold if the correct distance is greater.
+    /// This is for performace only
+    /// </summary>
+    /// <param name="a"></param>
+    /// <param name="b"></param>
+    /// <param name="threshold"></param>
+    /// <returns></returns>
+    double Get(T a, T b, double threshold);
+
+    /// <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/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/TSNEInterfaces/IHeap.cs
===================================================================
--- /branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/TSNEInterfaces/IHeap.cs	(revision 14414)
+++ /branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/TSNEInterfaces/IHeap.cs	(revision 14414)
@@ -0,0 +1,34 @@
+﻿#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;
+
+namespace HeuristicLab.Algorithms.DataAnalysis {
+  public interface IHeap<TK, TV> {
+    int Size { get; }
+    KeyValuePair<TK, TV> PeekMin();
+    TK PeekMinKey();
+    TV PeekMinValue();
+    void Insert(TK key, TV value);
+    void RemoveMin();
+  }
+}
Index: /branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/TSNEInterfaces/ISPTree.cs
===================================================================
--- /branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/TSNEInterfaces/ISPTree.cs	(revision 14414)
+++ /branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/TSNEInterfaces/ISPTree.cs	(revision 14414)
@@ -0,0 +1,39 @@
+﻿#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 ISPTree : IDeepCloneable {
+
+    double[,] Data { set; }
+    int GetDepth();
+    ISPTree GetParent();
+    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, double[] sumQ);
+    void ComputeEdgeForces(int[] rowP, int[] colP, double[] valP, int n, double[,] posF);
+
+  }
+}
Index: /branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/TSNEInterfaces/ITSNE.cs
===================================================================
--- /branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/TSNEInterfaces/ITSNE.cs	(revision 14414)
+++ /branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/TSNEInterfaces/ITSNE.cs	(revision 14414)
@@ -0,0 +1,28 @@
+﻿#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 ITSNE<in T> : IDeepCloneable {
+    double[,] Run(T[] data, int newDimensions, double perplexity, double theta);
+  }
+}
Index: /branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/TSNEInterfaces/IVPTree.cs
===================================================================
--- /branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/TSNEInterfaces/IVPTree.cs	(revision 14414)
+++ /branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/TSNEInterfaces/IVPTree.cs	(revision 14414)
@@ -0,0 +1,30 @@
+﻿#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 IVPTree<T> : IDeepCloneable {
+    void Create(IEnumerable<T> items);
+    void Search(T target, int k, out List<T> results, out List<double> distances);
+  }
+}
Index: /branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/TSNEUtils.cs
===================================================================
--- /branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/TSNEUtils.cs	(revision 14414)
+++ /branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/TSNEUtils.cs	(revision 14414)
@@ -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/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/VPTree.cs
===================================================================
--- /branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/VPTree.cs	(revision 14414)
+++ /branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/VPTree.cs	(revision 14414)
@@ -0,0 +1,223 @@
+﻿#region License Information
+/* HeuristicLab
+ * Copyright (C) 2002-2016 Heuristic and Evolutionary Algorithms Laboratory (HEAL)
+ *
+ * This file is part of HeuristicLab.
+ *
+ * HeuristicLab is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * HeuristicLab is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with HeuristicLab. If not, see <http://www.gnu.org/licenses/>.
+ */
+
+//Code is based on an implementation from Laurens van der Maaten
+
+/*
+*
+* Copyright (c) 2014, Laurens van der Maaten (Delft University of Technology)
+* All rights reserved.
+*
+* Redistribution and use in source and binary forms, with or without
+* modification, are permitted provided that the following conditions are met:
+* 1. Redistributions of source code must retain the above copyright
+*    notice, this list of conditions and the following disclaimer.
+* 2. Redistributions in binary form must reproduce the above copyright
+*    notice, this list of conditions and the following disclaimer in the
+*    documentation and/or other materials provided with the distribution.
+* 3. All advertising materials mentioning features or use of this software
+*    must display the following acknowledgement:
+*    This product includes software developed by the Delft University of Technology.
+* 4. Neither the name of the Delft University of Technology nor the names of
+*    its contributors may be used to endorse or promote products derived from
+*    this software without specific prior written permission.
+*
+* THIS SOFTWARE IS PROVIDED BY LAURENS VAN DER MAATEN ''AS IS'' AND ANY EXPRESS
+* OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
+* OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO
+* EVENT SHALL LAURENS VAN DER MAATEN BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+* SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+* PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR
+* BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+* CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING
+* IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY
+* OF SUCH DAMAGE.
+*
+*/
+#endregion
+
+using System;
+using System.Collections.Generic;
+using System.Linq;
+using HeuristicLab.Common;
+using HeuristicLab.Core;
+using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
+using HeuristicLab.Random;
+
+namespace HeuristicLab.Algorithms.DataAnalysis {
+  [StorableClass]
+  public class VPTree<T> : Item, IVPTree<T> where T : class, IDeepCloneable {
+    #region properties
+    [Storable]
+    private List<T> items;
+    [Storable]
+    private double tau;
+    [Storable]
+    private Node root;
+    [Storable]
+    private IDistance<T> distance;
+    #endregion
+
+    #region HLConstructors & Cloning
+    [StorableConstructor]
+    protected VPTree(bool deserializing) : base(deserializing) { }
+    protected VPTree(VPTree<T> original, Cloner cloner)
+      : base(original, cloner) {
+      items = original.items.Select(cloner.Clone).ToList();
+      tau = original.tau;
+      root = cloner.Clone(original.root);
+      distance = cloner.Clone(distance);
+    }
+    public override IDeepCloneable Clone(Cloner cloner) { return new VPTree<T>(this, cloner); }
+    public VPTree(IDistance<T> distance) {
+      root = null;
+      this.distance = distance;
+    }
+    #endregion
+
+    public void Create(IEnumerable<T> items) {
+      this.items = items.Select(x => x).ToList();
+      root = BuildFromPoints(0, this.items.Count);
+    }
+
+    public void Search(T target, int k, out List<T> results, out List<double> distances) {
+      IHeap<double, HeapItem> heap = new PriorityQueue<double, HeapItem>(double.MaxValue, double.MinValue, k);
+      tau = double.MaxValue;
+      Search(root, target, k, heap);
+      results = new List<T>();
+      distances = new List<double>();
+      while (heap.Size > 0) {
+        results.Add(items[heap.PeekMinValue().Index]);
+        distances.Add(heap.PeekMinValue().Dist);
+        heap.RemoveMin();
+      }
+      results.Reverse();
+      distances.Reverse();
+    }
+
+    private Node BuildFromPoints(int lower, int upper) {
+      if (upper == lower)      // indicates that we're done here!
+        return null;
+
+      // Lower index is center of current node
+      var node = new Node { index = lower };
+      var r = new MersenneTwister(); //outer behaviour does not change with the random seed => no need to take the IRandom from the algorithm 
+      if (upper - lower <= 1) return node; // if we did not arrive at leaf yet
+
+      // Choose an arbitrary point and move it to the start
+      var i = (int)(r.NextDouble() / 1 * (upper - lower - 1)) + lower;
+      items.Swap(lower, i);
+
+      // Partition around the median distance
+      var median = (upper + lower) / 2;
+      items.NthElement(lower + 1, upper - 1, median, distance.GetDistanceComparer(items[lower]));
+
+      // Threshold of the new node will be the distance to the median
+      node.threshold = distance.Get(items[lower], items[median]);
+
+      // Recursively build tree
+      node.index = lower;
+      node.left = BuildFromPoints(lower + 1, median);
+      node.right = BuildFromPoints(median, upper);
+
+      // Return result
+      return node;
+    }
+
+    private void Search(Node node, T target, int k, IHeap<double, HeapItem> heap) {
+      if (node == null) return;
+      var dist = distance.Get(items[node.index], target);
+      if (dist < tau) {
+        if (heap.Size == k) heap.RemoveMin();
+        heap.Insert(-dist, new HeapItem(node.index, dist));//TODO check if minheap or maxheap schould be used here
+        if (heap.Size == k) tau = heap.PeekMinValue().Dist;
+      }
+      if (node.left == null && node.right == null) return;
+
+      if (dist < node.threshold) {
+        if (dist - tau <= node.threshold) Search(node.left, target, k, heap);   // if there can still be neighbors inside the ball, recursively search left child first
+        if (dist + tau >= node.threshold) Search(node.right, target, k, heap);  // if there can still be neighbors outside the ball, recursively search right child
+      } else {
+        if (dist + tau >= node.threshold) Search(node.right, target, k, heap);  // if there can still be neighbors outside the ball, recursively search right child first
+        if (dist - tau <= node.threshold) Search(node.left, target, k, heap);   // if there can still be neighbors inside the ball, recursively search left child
+      }
+
+    }
+
+    [StorableClass]
+    private sealed class Node : Item {
+      [Storable]
+      public int index;
+      [Storable]
+      public double threshold;
+      [Storable]
+      public Node left;
+      [Storable]
+      public Node right;
+
+      #region HLConstructors & Cloning
+      [StorableConstructor]
+      private Node(bool deserializing) : base(deserializing) { }
+      private Node(Node original, Cloner cloner) : base(original, cloner) {
+        index = original.index;
+        threshold = original.threshold;
+        left = (Node)original.left.Clone(cloner);
+        right = (Node)original.right.Clone(cloner);
+      }
+      internal Node() {
+        index = 0;
+        threshold = 0;
+        left = null;
+        right = null;
+      }
+      public override IDeepCloneable Clone(Cloner cloner) {
+        return new Node(this, cloner);
+      }
+      #endregion
+    }
+
+    [StorableClass]
+    private sealed class HeapItem : Item, IComparable<HeapItem> {
+      [Storable]
+      public int Index;
+      [Storable]
+      public double Dist;
+
+      #region HLConstructors & Cloning
+      [StorableConstructor]
+      private HeapItem(bool deserializing) : base(deserializing) { }
+      private HeapItem(HeapItem original, Cloner cloner)
+      : base(original, cloner) {
+        Index = original.Index;
+        Dist = original.Dist;
+      }
+      public override IDeepCloneable Clone(Cloner cloner) { return new HeapItem(this, cloner); }
+      public HeapItem(int index, double dist) {
+        Index = index;
+        Dist = dist;
+      }
+      #endregion
+
+      public int CompareTo(HeapItem other) {
+        return Dist.CompareTo(Dist);
+      }
+    }
+  }
+}
