Index: branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/TSNE.cs
===================================================================
--- branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/TSNE.cs	(revision 14778)
+++ branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/TSNE.cs	(revision 14782)
@@ -1,3 +1,3 @@
-﻿#region License Information
+#region License Information
 /* HeuristicLab
  * Copyright (C) 2002-2016 Heuristic and Evolutionary Algorithms Laboratory (HEAL)
@@ -18,628 +18,265 @@
  * 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.Drawing;
 using System.Linq;
+using System.Threading;
 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>
+  /// t-distributed stochastic neighbourhood embedding (TSNE) projects the data in a low dimensional 
+  /// space to allow visual cluster identification.
+  /// </summary>
+  [Item("TSNE", "t-distributed stochastic neighbourhood embedding projects the data in a low " +
+                "dimensional space to allow visual cluster identification.")]
+  [Creatable(CreatableAttribute.Categories.DataAnalysis, Priority = 100)]
   [StorableClass]
-  public class TSNE<T> : DeepCloneable where T : class, IDeepCloneable {
-
-    private const string IterationResultName = "Iteration";
-    private const string ErrorResultName = "Error";
-    private const string ErrorPlotResultName = "ErrorPlot";
-    private const string ScatterPlotResultName = "Scatterplot";
-    private const string DataResultName = "Projected Data";
-
-    #region Properties
+  public sealed class TSNEAnalysis : BasicAlgorithm {
+    public override bool SupportsPause {
+      get { return false; }
+    }
+    public override Type ProblemType {
+      get { return typeof(IDataAnalysisProblem); }
+    }
+    public new IDataAnalysisProblem Problem {
+      get { return (IDataAnalysisProblem)base.Problem; }
+      set { base.Problem = value; }
+    }
+
+    #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";
+    private const string ClassesParameterName = "ClassNames";
+    private const string NormalizationParameterName = "Normalization";
+    #endregion
+
+    #region Parameterproperties
+    public IFixedValueParameter<DoubleValue> PerplexityParameter {
+      get { return Parameters[PerplexityParameterName] as IFixedValueParameter<DoubleValue>; }
+    }
+    public OptionalValueParameter<DoubleValue> ThetaParameter {
+      get { return Parameters[ThetaParameterName] as OptionalValueParameter<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>; }
+    }
+    public IFixedValueParameter<StringValue> ClassesParameter {
+      get { return Parameters[ClassesParameterName] as IFixedValueParameter<StringValue>; }
+    }
+    public IFixedValueParameter<BoolValue> NormalizationParameter {
+      get { return Parameters[NormalizationParameterName] as IFixedValueParameter<BoolValue>; }
+    }
+    #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 == null ? 0 : 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 == null ? 0 : EtaParameter.Value.Value;
+      }
+    }
+    public bool SetSeedRandomly {
+      get { return SetSeedRandomlyParameter.Value.Value; }
+    }
+    public uint Seed {
+      get { return (uint)SeedParameter.Value.Value; }
+    }
+    public string Classes {
+      get { return ClassesParameter.Value.Value; }
+    }
+    public bool Normalization {
+      get { return NormalizationParameter.Value.Value; }
+    }
     [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;
-    [Storable]
-    private Dictionary<string, List<int>> dataRowLookup;
-    [Storable]
-    private Dictionary<string, ScatterPlotDataRow> dataRows;
-    #endregion
-
-    #region Stopping 
-    public volatile bool Running;
-    #endregion
-
-    #region HLConstructors & Cloning
+    public TSNE<RealVector> tsne;
+    #endregion
+
+    #region Constructors & Cloning
     [StorableConstructor]
-    protected TSNE(bool 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);
-      dataRowLookup = original.dataRowLookup.ToDictionary(entry => entry.Key, entry => entry.Value.Select(x => x).ToList());
-      dataRows = original.dataRows.ToDictionary(entry => entry.Key, entry => cloner.Clone(entry.Value));
-    }
-    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, Dictionary<string, List<int>> dataRowLookup = null, Dictionary<string, ScatterPlotDataRow> dataRows = null) {
-      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;
-      this.dataRowLookup = dataRowLookup;
-      if (dataRows != null) this.dataRows = dataRows;
-      else { this.dataRows = new Dictionary<string, ScatterPlotDataRow>(); }
-    }
-    #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!");
-      SetUpResults(data);
-      Running = true;
-      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;
-      double[,] p = null;
-      int[] rowP = null;
-      int[] colP = null;
-      double[] valP = null;
-      var rand = new NormalDistributedRandom(random, 0, 1);
-
-      //Calculate Similarities
-      if (exact) p = CalculateExactSimilarites(data, perplexity);
-      else CalculateApproximateSimilarities(data, perplexity, out rowP, out colP, out valP);
-
-      // 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;
-
-      // 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 && Running; iter++) {
-        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;
-
-        Analyze(exact, iter, p, rowP, colP, valP, newData, noDatapoints, newDimensions, theta);
-      }
-      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 SetUpResults(IReadOnlyCollection<T> data) {
-      if (dataRowLookup == null) dataRowLookup = new Dictionary<string, List<int>> { { "Data", Enumerable.Range(0, data.Count).ToList() } };
-      if (results == null) return;
-
-      if (!results.ContainsKey(IterationResultName)) results.Add(new Result(IterationResultName, new IntValue(0)));
-      else ((IntValue)results[IterationResultName].Value).Value = 0;
-
-      if (!results.ContainsKey(ErrorResultName)) results.Add(new Result(ErrorResultName, new DoubleValue(0)));
-      else ((DoubleValue)results[ErrorResultName].Value).Value = 0;
-
-      if (!results.ContainsKey(ErrorPlotResultName)) results.Add(new Result(ErrorPlotResultName, new DataTable(ErrorPlotResultName, "Development of errors during Gradiant descent")));
-      else results[ErrorPlotResultName].Value = new DataTable(ErrorPlotResultName, "Development of errors during Gradiant descent");
-
-      var plot = results[ErrorPlotResultName].Value as DataTable;
-      if (plot == null) throw new ArgumentException("could not create/access Error-DataTable in Results-Collection");
-
-      if (!plot.Rows.ContainsKey("errors")) plot.Rows.Add(new DataRow("errors"));
-      plot.Rows["errors"].Values.Clear();
-
-      results.Add(new Result(ScatterPlotResultName, "Plot of the projected data", new ScatterPlot(DataResultName, "")));
-      results.Add(new Result(DataResultName, "Projected Data", new DoubleMatrix()));
-
-    }
-    private void Analyze(bool exact, int iter, double[,] p, int[] rowP, int[] colP, double[] valP, double[,] newData, int noDatapoints, int newDimensions, double theta) {
-      if (results == null) return;
-      var plot = results[ErrorPlotResultName].Value as DataTable;
-      if (plot == null) throw new ArgumentException("Could not create/access error data table in results collection. Was it removed by some effect?");
-      var errors = plot.Rows["errors"].Values;
-      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();
-
-      var ndata = Normalize(newData);
-      results[DataResultName].Value = new DoubleMatrix(ndata);
-      var splot = results[ScatterPlotResultName].Value as ScatterPlot;
-      FillScatterPlot(ndata, splot);
-
-
-    }
-    private void FillScatterPlot(double[,] lowDimData, ScatterPlot plot) {
-      foreach (var rowName in dataRowLookup.Keys) {
-        if (!plot.Rows.ContainsKey(rowName))
-          plot.Rows.Add(dataRows.ContainsKey(rowName) ? dataRows[rowName] : new ScatterPlotDataRow(rowName, "", new List<Point2D<double>>()));
-        plot.Rows[rowName].Points.Replace(dataRowLookup[rowName].Select(i => new Point2D<double>(lowDimData[i, 0], lowDimData[i, 1])));
-      }
-    }
-    private static double[,] Normalize(double[,] data) {
-      var max = new double[data.GetLength(1)];
-      var min = new double[data.GetLength(1)];
-      var res = new double[data.GetLength(0), data.GetLength(1)];
-      for (var i = 0; i < max.Length; i++) max[i] = min[i] = data[0, i];
-      for (var i = 0; i < data.GetLength(0); i++)
-        for (var j = 0; j < data.GetLength(1); j++) {
-          var v = data[i, j];
-          max[j] = Math.Max(max[j], v);
-          min[j] = Math.Min(min[j], v);
-        }
-      for (var i = 0; i < data.GetLength(0); i++) {
-        for (var j = 0; j < data.GetLength(1); j++) {
-          res[i, j] = (data[i, j] - (max[j] + min[j]) / 2) / (max[j] - min[j]);
-        }
-      }
-      return res;
-    }
-    private void CalculateApproximateSimilarities(T[] data, double perplexity, out int[] rowP, out int[] colP, out double[] valP) {
-      // Compute asymmetric pairwise input similarities
-      ComputeGaussianPerplexity(data, data.Length, out rowP, out colP, out valP, perplexity, (int)(3 * perplexity));
-      // Symmetrize input similarities
-      int[] sRowP, symColP;
-      double[] sValP;
-      SymmetrizeMatrix(rowP, colP, valP, out sRowP, out symColP, out sValP);
-      rowP = sRowP;
-      colP = symColP;
-      valP = sValP;
-      var sumP = .0;
-      for (var i = 0; i < rowP[data.Length]; i++) sumP += valP[i];
-      for (var i = 0; i < rowP[data.Length]; i++) valP[i] /= sumP;
-    }
-    private double[,] CalculateExactSimilarites(T[] data, double perplexity) {
-      // Compute similarities
-      var p = new double[data.Length, data.Length];
-      ComputeGaussianPerplexity(data, data.Length, p, perplexity);
-      // Symmetrize input similarities
-      for (var n = 0; n < data.Length; n++) {
-        for (var m = n + 1; m < data.Length; m++) {
-          p[n, m] += p[m, n];
-          p[m, n] = p[n, m];
-        }
-      }
-      var sumP = .0;
-      for (var i = 0; i < data.Length; i++) for (var j = 0; j < data.Length; j++) sumP += p[i, j];
-      for (var i = 0; i < data.Length; i++) for (var j = 0; j < data.Length; j++) p[i, j] /= sumP;
-      return p;
-    }
-
-    private 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;
-            }
+    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. Recommended Value is Floor(number of points /3) or lower", new DoubleValue(25)));
+      Parameters.Add(new OptionalValueParameter<DoubleValue>(ThetaParameterName, "Value describing how much appoximated gradients my differ from exact gradients. Set to 0 for exact calculation and in [0,1] otherwise \n CAUTION: exact calculation of forces requires building a non-sparse N*N matrix where N is the number of data points\n This may exceed memory limitations", 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(0)));
+      Parameters.Add(new FixedValueParameter<IntValue>(MomentumSwitchIterationParameterName, "Number of iterations after which the momentum in the gradient descent is switched", new IntValue(0)));
+      Parameters.Add(new FixedValueParameter<DoubleValue>(InitialMomentumParameterName, "The initial momentum in the gradient descent", new DoubleValue(0.5)));
+      Parameters.Add(new FixedValueParameter<DoubleValue>(FinalMomentumParameterName, "The final momentum", new DoubleValue(0.8)));
+      Parameters.Add(new FixedValueParameter<DoubleValue>(EtaParameterName, "Gradient Descent learning rate", new DoubleValue(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)));
+      Parameters.Add(new FixedValueParameter<StringValue>(ClassesParameterName, "name of the column specifying the class lables of each data point. \n if the lable column can not be found Training/Test is used as labels", new StringValue("none")));
+      Parameters.Add(new FixedValueParameter<BoolValue>(NormalizationParameterName, "Wether the data should be zero centered and have variance of 1 for each variable, so different scalings are ignored", new BoolValue(true)));
+
+      MomentumSwitchIterationParameter.Hidden = true;
+      InitialMomentumParameter.Hidden = true;
+      FinalMomentumParameter.Hidden = true;
+      StopLyingIterationParameter.Hidden = true;
+      EtaParameter.Hidden = true;
+    }
+    #endregion
+
+    public override void Stop() {
+      base.Stop();
+      if(tsne != null) tsne.Running = false;
+    }
+
+    protected override void Run(CancellationToken cancellationToken) {
+      var dataRowNames = new Dictionary<string, List<int>>();
+      var rows = new Dictionary<string, ScatterPlotDataRow>();
+      var problemData = Problem.ProblemData;
+
+      //color datapoints acording to Classes-Variable (be it double or string)
+      if(problemData.Dataset.VariableNames.Contains(Classes)) {
+        if((problemData.Dataset as Dataset).VariableHasType<string>(Classes)) {
+          var classes = problemData.Dataset.GetStringValues(Classes).ToArray();
+          for(var i = 0; i < classes.Length; i++) {
+            if(!dataRowNames.ContainsKey(classes[i])) dataRowNames.Add(classes[i], new List<int>());
+            dataRowNames[classes[i]].Add(i);
           }
-
-          // 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;
-            }
+        } else if((problemData.Dataset as Dataset).VariableHasType<double>(Classes)) {
+          var classValues = problemData.Dataset.GetDoubleValues(Classes).ToArray();
+          var max = classValues.Max() + 0.1;
+          var min = classValues.Min() - 0.1;
+          const int contours = 8;
+          for(var i = 0; i < contours; i++) {
+            var contourname = GetContourName(i, min, max, contours);
+            dataRowNames.Add(contourname, new List<int>());
+            rows.Add(contourname, new ScatterPlotDataRow(contourname, "", new List<Point2D<double>>()));
+            rows[contourname].VisualProperties.Color = GetHeatMapColor(i, contours);
+            rows[contourname].VisualProperties.PointSize = i + 3;
           }
-
-          // 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;
+          for(var i = 0; i < classValues.Length; i++) {
+            dataRowNames[GetContourName(classValues[i], min, max, contours)].Add(i);
           }
         }
-      }
-    }
-    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 (var 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));
-        }
+      } else {
+        dataRowNames.Add("Training", problemData.TrainingIndices.ToList());
+        dataRowNames.Add("Test", problemData.TestIndices.ToList());
+      }
+
+      //Set up and run TSNE
+      if(SetSeedRandomly) SeedParameter.Value.Value = new System.Random().Next();
+      var random = new MersenneTwister(Seed);
+      tsne = new TSNE<RealVector>(Distance, random, Results, MaxIterations, StopLyingIteration, MomentumSwitchIteration, InitialMomentum, FinalMomentum, Eta, dataRowNames, rows);
+      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());
+      if(Normalization) data = NormalizeData(data);
+      tsne.Run(data, NewDimensions, Perplexity, Theta);
+    }
+
+    private static RealVector[] NormalizeData(IReadOnlyList<RealVector> data) {
+      var n = data[0].Length;
+      var mean = new double[n];
+      var sd = new double[n];
+      var nData = new RealVector[data.Count];
+      for(var i = 0; i < n; i++) {
+        var i1 = i;
+        sd[i] = Enumerable.Range(0, data.Count).Select(x => data[x][i1]).StandardDeviation();
+        mean[i] = Enumerable.Range(0, data.Count).Select(x => data[x][i1]).Average();
+      }
+      for(var i = 0; i < data.Count; i++) {
+        nData[i] = new RealVector(n);
+        for(var j = 0; j < n; j++) nData[i][j] = (data[i][j] - mean[j]) / sd[j];
+      }
+      return nData;
+    }
+    private static Color GetHeatMapColor(int contourNr, int noContours) {
+      var q = (double)contourNr / noContours;  // q in [0,1]
+      var c = q < 0.5 ? Color.FromArgb((int)(q * 2 * 255), 255, 0) : Color.FromArgb(255, (int)((1 - q) * 2 * 255), 0);
       return c;
     }
-    private static 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
+    private static string GetContourName(double value, double min, double max, int noContours) {
+      var size = (max - min) / noContours;
+      var contourNr = (int)((value - min) / size);
+      return GetContourName(contourNr, min, max, noContours);
+    }
+    private static string GetContourName(int i, double min, double max, int noContours) {
+      var size = (max - min) / noContours;
+      return "[" + (min + i * size) + ";" + (min + (i + 1) * size) + ")";
+    }
   }
 }
Index: branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/TSNEAnalysis.cs
===================================================================
--- branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/TSNEAnalysis.cs	(revision 14778)
+++ 	(revision )
@@ -1,282 +1,0 @@
-#region License Information
-/* HeuristicLab
- * Copyright (C) 2002-2016 Heuristic and Evolutionary Algorithms Laboratory (HEAL)
- *
- * This file is part of HeuristicLab.
- *
- * HeuristicLab is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * HeuristicLab is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with HeuristicLab. If not, see <http://www.gnu.org/licenses/>.
- */
-#endregion
-
-using System;
-using System.Collections.Generic;
-using System.Drawing;
-using System.Linq;
-using System.Threading;
-using HeuristicLab.Analysis;
-using HeuristicLab.Common;
-using HeuristicLab.Core;
-using HeuristicLab.Data;
-using HeuristicLab.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>
-  /// t-distributed stochastic neighbourhood embedding (TSNE) projects the data in a low dimensional 
-  /// space to allow visual cluster identification.
-  /// </summary>
-  [Item("TSNE", "t-distributed stochastic neighbourhood embedding projects the data in a low " +
-                "dimensional space to allow visual cluster identification.")]
-  [Creatable(CreatableAttribute.Categories.DataAnalysis, Priority = 100)]
-  [StorableClass]
-  public sealed class TSNEAnalysis : BasicAlgorithm {
-    public override bool SupportsPause {
-      get { return false; }
-    }
-    public override Type ProblemType {
-      get { return typeof(IDataAnalysisProblem); }
-    }
-    public new IDataAnalysisProblem Problem {
-      get { return (IDataAnalysisProblem)base.Problem; }
-      set { base.Problem = value; }
-    }
-
-    #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";
-    private const string ClassesParameterName = "ClassNames";
-    private const string NormalizationParameterName = "Normalization";
-    #endregion
-
-    #region Parameterproperties
-    public IFixedValueParameter<DoubleValue> PerplexityParameter {
-      get { return Parameters[PerplexityParameterName] as IFixedValueParameter<DoubleValue>; }
-    }
-    public OptionalValueParameter<DoubleValue> ThetaParameter {
-      get { return Parameters[ThetaParameterName] as OptionalValueParameter<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>; }
-    }
-    public IFixedValueParameter<StringValue> ClassesParameter {
-      get { return Parameters[ClassesParameterName] as IFixedValueParameter<StringValue>; }
-    }
-    public IFixedValueParameter<BoolValue> NormalizationParameter {
-      get { return Parameters[NormalizationParameterName] as IFixedValueParameter<BoolValue>; }
-    }
-    #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 == null ? 0 : 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 == null ? 0 : EtaParameter.Value.Value;
-      }
-    }
-    public bool SetSeedRandomly {
-      get { return SetSeedRandomlyParameter.Value.Value; }
-    }
-    public uint Seed {
-      get { return (uint)SeedParameter.Value.Value; }
-    }
-    public string Classes {
-      get { return ClassesParameter.Value.Value; }
-    }
-    public bool Normalization {
-      get { return NormalizationParameter.Value.Value; }
-    }
-    [Storable]
-    public TSNE<RealVector> tsne;
-    #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. Recommended Value is Floor(number of points /3) or lower", new DoubleValue(25)));
-      Parameters.Add(new OptionalValueParameter<DoubleValue>(ThetaParameterName, "Value describing how much appoximated gradients my differ from exact gradients. Set to 0 for exact calculation and in [0,1] otherwise \n CAUTION: exact calculation of forces requires building a non-sparse N*N matrix where N is the number of data points\n This may exceed memory limitations", 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(0)));
-      Parameters.Add(new FixedValueParameter<IntValue>(MomentumSwitchIterationParameterName, "Number of iterations after which the momentum in the gradient descent is switched", new IntValue(0)));
-      Parameters.Add(new FixedValueParameter<DoubleValue>(InitialMomentumParameterName, "The initial momentum in the gradient descent", new DoubleValue(0.5)));
-      Parameters.Add(new FixedValueParameter<DoubleValue>(FinalMomentumParameterName, "The final momentum", new DoubleValue(0.8)));
-      Parameters.Add(new FixedValueParameter<DoubleValue>(EtaParameterName, "Gradient Descent learning rate", new DoubleValue(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)));
-      Parameters.Add(new FixedValueParameter<StringValue>(ClassesParameterName, "name of the column specifying the class lables of each data point. \n if the lable column can not be found Training/Test is used as labels", new StringValue("none")));
-      Parameters.Add(new FixedValueParameter<BoolValue>(NormalizationParameterName, "Wether the data should be zero centered and have variance of 1 for each variable, so different scalings are ignored", new BoolValue(true)));
-
-      MomentumSwitchIterationParameter.Hidden = true;
-      InitialMomentumParameter.Hidden = true;
-      FinalMomentumParameter.Hidden = true;
-      StopLyingIterationParameter.Hidden = true;
-      EtaParameter.Hidden = true;
-    }
-    #endregion
-
-    public override void Stop() {
-      base.Stop();
-      if(tsne != null) tsne.Running = false;
-    }
-
-    protected override void Run(CancellationToken cancellationToken) {
-      var dataRowNames = new Dictionary<string, List<int>>();
-      var rows = new Dictionary<string, ScatterPlotDataRow>();
-      var problemData = Problem.ProblemData;
-
-      //color datapoints acording to Classes-Variable (be it double or string)
-      if(problemData.Dataset.VariableNames.Contains(Classes)) {
-        if((problemData.Dataset as Dataset).VariableHasType<string>(Classes)) {
-          var classes = problemData.Dataset.GetStringValues(Classes).ToArray();
-          for(var i = 0; i < classes.Length; i++) {
-            if(!dataRowNames.ContainsKey(classes[i])) dataRowNames.Add(classes[i], new List<int>());
-            dataRowNames[classes[i]].Add(i);
-          }
-        } else if((problemData.Dataset as Dataset).VariableHasType<double>(Classes)) {
-          var classValues = problemData.Dataset.GetDoubleValues(Classes).ToArray();
-          var max = classValues.Max() + 0.1;
-          var min = classValues.Min() - 0.1;
-          const int contours = 8;
-          for(var i = 0; i < contours; i++) {
-            var contourname = GetContourName(i, min, max, contours);
-            dataRowNames.Add(contourname, new List<int>());
-            rows.Add(contourname, new ScatterPlotDataRow(contourname, "", new List<Point2D<double>>()));
-            rows[contourname].VisualProperties.Color = GetHeatMapColor(i, contours);
-            rows[contourname].VisualProperties.PointSize = i + 3;
-          }
-          for(var i = 0; i < classValues.Length; i++) {
-            dataRowNames[GetContourName(classValues[i], min, max, contours)].Add(i);
-          }
-        }
-      } else {
-        dataRowNames.Add("Training", problemData.TrainingIndices.ToList());
-        dataRowNames.Add("Test", problemData.TestIndices.ToList());
-      }
-
-      //Set up and run TSNE
-      if(SetSeedRandomly) SeedParameter.Value.Value = new System.Random().Next();
-      var random = new MersenneTwister(Seed);
-      tsne = new TSNE<RealVector>(Distance, random, Results, MaxIterations, StopLyingIteration, MomentumSwitchIteration, InitialMomentum, FinalMomentum, Eta, dataRowNames, rows);
-      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());
-      if(Normalization) data = NormalizeData(data);
-      tsne.Run(data, NewDimensions, Perplexity, Theta);
-    }
-
-    private static RealVector[] NormalizeData(IReadOnlyList<RealVector> data) {
-      var n = data[0].Length;
-      var mean = new double[n];
-      var sd = new double[n];
-      var nData = new RealVector[data.Count];
-      for(var i = 0; i < n; i++) {
-        var i1 = i;
-        sd[i] = Enumerable.Range(0, data.Count).Select(x => data[x][i1]).StandardDeviation();
-        mean[i] = Enumerable.Range(0, data.Count).Select(x => data[x][i1]).Average();
-      }
-      for(var i = 0; i < data.Count; i++) {
-        nData[i] = new RealVector(n);
-        for(var j = 0; j < n; j++) nData[i][j] = (data[i][j] - mean[j]) / sd[j];
-      }
-      return nData;
-    }
-    private static Color GetHeatMapColor(int contourNr, int noContours) {
-      var q = (double)contourNr / noContours;  // q in [0,1]
-      var c = q < 0.5 ? Color.FromArgb((int)(q * 2 * 255), 255, 0) : Color.FromArgb(255, (int)((1 - q) * 2 * 255), 0);
-      return c;
-    }
-    private static string GetContourName(double value, double min, double max, int noContours) {
-      var size = (max - min) / noContours;
-      var contourNr = (int)((value - min) / size);
-      return GetContourName(contourNr, min, max, noContours);
-    }
-    private static string GetContourName(int i, double min, double max, int noContours) {
-      var size = (max - min) / noContours;
-      return "[" + (min + i * size) + ";" + (min + (i + 1) * size) + ")";
-    }
-  }
-}
Index: branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/TSNEStatic.cs
===================================================================
--- branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/TSNEStatic.cs	(revision 14782)
+++ branches/TSNE/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/TSNEStatic.cs	(revision 14782)
@@ -0,0 +1,645 @@
+﻿#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.Analysis;
+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> : DeepCloneable where T : class, IDeepCloneable {
+
+    private const string IterationResultName = "Iteration";
+    private const string ErrorResultName = "Error";
+    private const string ErrorPlotResultName = "ErrorPlot";
+    private const string ScatterPlotResultName = "Scatterplot";
+    private const string DataResultName = "Projected Data";
+
+    #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;
+    [Storable]
+    private Dictionary<string, List<int>> dataRowLookup;
+    [Storable]
+    private Dictionary<string, ScatterPlotDataRow> dataRows;
+    #endregion
+
+    #region Stopping 
+    public volatile bool Running;
+    #endregion
+
+    #region HLConstructors & Cloning
+    [StorableConstructor]
+    protected TSNE(bool 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);
+      dataRowLookup = original.dataRowLookup.ToDictionary(entry => entry.Key, entry => entry.Value.Select(x => x).ToList());
+      dataRows = original.dataRows.ToDictionary(entry => entry.Key, entry => cloner.Clone(entry.Value));
+    }
+    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, Dictionary<string, List<int>> dataRowLookup = null, Dictionary<string, ScatterPlotDataRow> dataRows = null) {
+      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;
+      this.dataRowLookup = dataRowLookup;
+      if (dataRows != null) this.dataRows = dataRows;
+      else { this.dataRows = new Dictionary<string, ScatterPlotDataRow>(); }
+    }
+    #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!");
+      SetUpResults(data);
+      Running = true;
+      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;
+      double[,] p = null;
+      int[] rowP = null;
+      int[] colP = null;
+      double[] valP = null;
+      var rand = new NormalDistributedRandom(random, 0, 1);
+
+      //Calculate Similarities
+      if (exact) p = CalculateExactSimilarites(data, perplexity);
+      else CalculateApproximateSimilarities(data, perplexity, out rowP, out colP, out valP);
+
+      // 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;
+
+      // 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 && Running; iter++) {
+        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;
+
+        Analyze(exact, iter, p, rowP, colP, valP, newData, noDatapoints, newDimensions, theta);
+      }
+      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 SetUpResults(IReadOnlyCollection<T> data) {
+      if (dataRowLookup == null) dataRowLookup = new Dictionary<string, List<int>> { { "Data", Enumerable.Range(0, data.Count).ToList() } };
+      if (results == null) return;
+
+      if (!results.ContainsKey(IterationResultName)) results.Add(new Result(IterationResultName, new IntValue(0)));
+      else ((IntValue)results[IterationResultName].Value).Value = 0;
+
+      if (!results.ContainsKey(ErrorResultName)) results.Add(new Result(ErrorResultName, new DoubleValue(0)));
+      else ((DoubleValue)results[ErrorResultName].Value).Value = 0;
+
+      if (!results.ContainsKey(ErrorPlotResultName)) results.Add(new Result(ErrorPlotResultName, new DataTable(ErrorPlotResultName, "Development of errors during Gradiant descent")));
+      else results[ErrorPlotResultName].Value = new DataTable(ErrorPlotResultName, "Development of errors during Gradiant descent");
+
+      var plot = results[ErrorPlotResultName].Value as DataTable;
+      if (plot == null) throw new ArgumentException("could not create/access Error-DataTable in Results-Collection");
+
+      if (!plot.Rows.ContainsKey("errors")) plot.Rows.Add(new DataRow("errors"));
+      plot.Rows["errors"].Values.Clear();
+
+      results.Add(new Result(ScatterPlotResultName, "Plot of the projected data", new ScatterPlot(DataResultName, "")));
+      results.Add(new Result(DataResultName, "Projected Data", new DoubleMatrix()));
+
+    }
+    private void Analyze(bool exact, int iter, double[,] p, int[] rowP, int[] colP, double[] valP, double[,] newData, int noDatapoints, int newDimensions, double theta) {
+      if (results == null) return;
+      var plot = results[ErrorPlotResultName].Value as DataTable;
+      if (plot == null) throw new ArgumentException("Could not create/access error data table in results collection. Was it removed by some effect?");
+      var errors = plot.Rows["errors"].Values;
+      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();
+
+      var ndata = Normalize(newData);
+      results[DataResultName].Value = new DoubleMatrix(ndata);
+      var splot = results[ScatterPlotResultName].Value as ScatterPlot;
+      FillScatterPlot(ndata, splot);
+
+
+    }
+    private void FillScatterPlot(double[,] lowDimData, ScatterPlot plot) {
+      foreach (var rowName in dataRowLookup.Keys) {
+        if (!plot.Rows.ContainsKey(rowName))
+          plot.Rows.Add(dataRows.ContainsKey(rowName) ? dataRows[rowName] : new ScatterPlotDataRow(rowName, "", new List<Point2D<double>>()));
+        plot.Rows[rowName].Points.Replace(dataRowLookup[rowName].Select(i => new Point2D<double>(lowDimData[i, 0], lowDimData[i, 1])));
+      }
+    }
+    private static double[,] Normalize(double[,] data) {
+      var max = new double[data.GetLength(1)];
+      var min = new double[data.GetLength(1)];
+      var res = new double[data.GetLength(0), data.GetLength(1)];
+      for (var i = 0; i < max.Length; i++) max[i] = min[i] = data[0, i];
+      for (var i = 0; i < data.GetLength(0); i++)
+        for (var j = 0; j < data.GetLength(1); j++) {
+          var v = data[i, j];
+          max[j] = Math.Max(max[j], v);
+          min[j] = Math.Min(min[j], v);
+        }
+      for (var i = 0; i < data.GetLength(0); i++) {
+        for (var j = 0; j < data.GetLength(1); j++) {
+          res[i, j] = (data[i, j] - (max[j] + min[j]) / 2) / (max[j] - min[j]);
+        }
+      }
+      return res;
+    }
+    private void CalculateApproximateSimilarities(T[] data, double perplexity, out int[] rowP, out int[] colP, out double[] valP) {
+      // Compute asymmetric pairwise input similarities
+      ComputeGaussianPerplexity(data, data.Length, out rowP, out colP, out valP, perplexity, (int)(3 * perplexity));
+      // Symmetrize input similarities
+      int[] sRowP, symColP;
+      double[] sValP;
+      SymmetrizeMatrix(rowP, colP, valP, out sRowP, out symColP, out sValP);
+      rowP = sRowP;
+      colP = symColP;
+      valP = sValP;
+      var sumP = .0;
+      for (var i = 0; i < rowP[data.Length]; i++) sumP += valP[i];
+      for (var i = 0; i < rowP[data.Length]; i++) valP[i] /= sumP;
+    }
+    private double[,] CalculateExactSimilarites(T[] data, double perplexity) {
+      // Compute similarities
+      var p = new double[data.Length, data.Length];
+      ComputeGaussianPerplexity(data, data.Length, p, perplexity);
+      // Symmetrize input similarities
+      for (var n = 0; n < data.Length; n++) {
+        for (var m = n + 1; m < data.Length; m++) {
+          p[n, m] += p[m, n];
+          p[m, n] = p[n, m];
+        }
+      }
+      var sumP = .0;
+      for (var i = 0; i < data.Length; i++) for (var j = 0; j < data.Length; j++) sumP += p[i, j];
+      for (var i = 0; i < data.Length; i++) for (var j = 0; j < data.Length; j++) p[i, j] /= sumP;
+      return p;
+    }
+
+    private 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 (var 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
+  }
+}
