Free cookie consent management tool by TermsFeed Policy Generator

Changeset 15207


Ignore:
Timestamp:
07/12/17 15:17:41 (7 years ago)
Author:
bwerth
Message:

#2700 worked on TSNE (mostly removing unused code)

Location:
trunk/sources/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE
Files:
1 added
1 deleted
8 edited
1 moved

Legend:

Unmodified
Added
Removed
  • trunk/sources/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/Distances/CosineDistance.cs

    r15206 r15207  
    3434  /// </summary>
    3535  [StorableClass]
    36   [Item("InnerProductDistance", "The angluar distance as defined as a normalized distance measure dependent on the angle between two vectors.\nIt is designed for vectors with all positive coordinates")]
    37   public class InnerProductDistance : DistanceBase<IEnumerable<double>> {
     36  [Item("CosineDistance", "The angluar distance as defined as a normalized distance measure dependent on the angle between two vectors.\nIt is designed for vectors with all positive coordinates")]
     37  public class CosineDistance : DistanceBase<IEnumerable<double>> {
    3838
    39     #region HLConstructors
     39    #region HLConstructors & Cloning
    4040    [StorableConstructor]
    41     protected InnerProductDistance(bool deserializing) : base(deserializing) { }
    42     protected InnerProductDistance(InnerProductDistance original, Cloner cloner)
     41    protected CosineDistance(bool deserializing) : base(deserializing) { }
     42    protected CosineDistance(CosineDistance original, Cloner cloner)
    4343      : base(original, cloner) { }
    44     public InnerProductDistance() { }
     44    public CosineDistance() { }
    4545    public override IDeepCloneable Clone(Cloner cloner) {
    46       return new InnerProductDistance(this, cloner);
     46      return new CosineDistance(this, cloner);
    4747    }
    4848    #endregion
    4949
    5050    #region statics
    51     public static double GetDistance(IEnumerable<double> point1, IEnumerable<double> point2) {
    52       var xs = point1.GetEnumerator();
    53       var ys = point2.GetEnumerator();
    54       var sum = 0.0;
    55       while(xs.MoveNext() & ys.MoveNext()) {
    56         if(xs.Current < 0 || ys.Current < 0) throw new ArgumentException("Inner product distance is only defined for vectors with non-negative elements");
    57         sum += xs.Current * ys.Current;
     51    public static double GetDistance(IReadOnlyList<double> point1, IReadOnlyList<double> point2) {
     52      if (point1.Count != point2.Count) throw new ArgumentException("Cosine distance not defined on vectors of different length");
     53      var innerprod = 0.0;
     54      var length1 = 0.0;
     55      var length2 = 0.0;
     56
     57      for (var i = 0; i < point1.Count; i++) {
     58        double d1 = point1[i], d2 = point2[i];
     59        innerprod += d1 * d2;
     60        length1 += d1 * d1;
     61        length2 += d2 * d2;
    5862      }
    59       if(xs.MoveNext() || ys.MoveNext()) throw new ArgumentException("Enumerables contain a different number of elements");
    60       return sum;
     63      var l = Math.Sqrt(length1 * length2);
     64      if (l.IsAlmost(0)) throw new ArgumentException("Cosine distance is not defined on vectors of length 0");
     65      return 1 - innerprod / l;
    6166    }
    6267    #endregion
    6368    public override double Get(IEnumerable<double> a, IEnumerable<double> b) {
    64       return GetDistance(a, b);
     69      return GetDistance(a.ToArray(), b.ToArray());
    6570    }
    6671  }
  • trunk/sources/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/Distances/DistanceBase.cs

    r14911 r15207  
    3030  public abstract class DistanceBase<T> : Item, IDistance<T> {
    3131
    32     #region HLConstructors & Boilerplate
     32    #region HLConstructors & Cloning
    3333    [StorableConstructor]
    3434    protected DistanceBase(bool deserializing) : base(deserializing) { }
     
    4242      return new DistanceComparer(item, this);
    4343    }
    44 
    4544
    4645    public double Get(object x, object y) {
  • trunk/sources/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/Distances/EuclideanDistance.cs

    r14784 r15207  
    3232  public class EuclideanDistance : DistanceBase<IEnumerable<double>> {
    3333
    34     #region HLConstructors & Boilerplate
     34    #region HLConstructors & Cloning
    3535    [StorableConstructor]
    3636    protected EuclideanDistance(bool deserializing) : base(deserializing) { }
     
    4040    #endregion
    4141
    42     public static double GetDistance(double[] point1, double[] point2) {
    43       if (point1.Length != point2.Length) throw new ArgumentException("Euclidean distance not defined on vectors of different length");
    44       return Math.Sqrt(point1.Zip(point2, (a1, b1) => (a1 - b1) * (a1 - b1)).Sum());
     42    public static double GetDistance(IReadOnlyList<double> point1, IReadOnlyList<double> point2) {
     43      if (point1.Count != point2.Count) throw new ArgumentException("Euclidean distance not defined on vectors of different length");
     44      var sum = 0.0;
     45      for (var i = 0; i < point1.Count; i++) {
     46        var d = point1[i] - point2[i];
     47        sum += d * d;
     48      }
     49
     50      return Math.Sqrt(sum);
    4551    }
    4652
  • trunk/sources/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/Distances/ManhattanDistance.cs

    r14911 r15207  
    3232  public class ManhattanDistance : DistanceBase<IEnumerable<double>> {
    3333
    34     #region HLConstructors & Boilerplate
     34    #region HLConstructors & Cloning
    3535    [StorableConstructor]
    3636    protected ManhattanDistance(bool deserializing) : base(deserializing) { }
     
    4747    public static double GetDistance(double[] point1, double[] point2) {
    4848      if (point1.Length != point2.Length) throw new ArgumentException("Manhattan distance not defined on vectors of different length");
    49       return point1.Zip(point2, (a1, b1) => Math.Abs(a1 - b1)).Sum();
     49      var sum = 0.0;
     50      for (var i = 0; i < point1.Length; i++)
     51        sum += Math.Abs(point1[i] + point2[i]);
     52      return sum;
    5053    }
    5154
  • trunk/sources/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/PriorityQueue.cs

    r14785 r15207  
    2525
    2626namespace HeuristicLab.Algorithms.DataAnalysis {
    27   // TODO: move to HeuristicLab.Collections
    2827
    2928  // Code Taken from branch: RoutePlanning, Heap4
     
    6362
    6463    public KeyValuePair<TK, TV> PeekMin() {
    65       if(size == 0) {
    66         throw new InvalidOperationException("Heap is empty");
    67       }
     64      if (size == 0) throw new InvalidOperationException("Heap is empty");
    6865      return new KeyValuePair<TK, TV>(data[1].Key, data[1].Value);
    6966    }
     
    8683      size = sz - 1;
    8784      finalLayerDist++;
    88       if(finalLayerDist == finalLayerSize) {
     85      if (finalLayerDist == finalLayerSize) {
    8986        finalLayerSize >>= 2;
    9087        finalLayerDist = 0;
    9188      }
    9289
    93       while(succ < sz) {
     90      while (succ < sz) {
    9491        var minKey = data[succ].Key;
    9592        var delta = 0;
    9693
    97         var otherKey = data[succ + 1].Key;
    98         if(otherKey.CompareTo(minKey) < 0) {
     94        for (var i = 1; i <= 3; i++) {
     95          var otherKey = data[succ + i].Key;
     96          if (otherKey.CompareTo(minKey) >= 0) continue;
    9997          minKey = otherKey;
    100           delta = 1;
     98          delta = i;
    10199        }
    102         otherKey = data[succ + 2].Key;
    103         if(otherKey.CompareTo(minKey) < 0) {
    104           minKey = otherKey;
    105           delta = 2;
    106         }
    107         otherKey = data[succ + 3].Key;
    108         if(otherKey.CompareTo(minKey) < 0) {
    109           minKey = otherKey;
    110           delta = 3;
    111         }
     100       
    112101        succ += delta;
    113102        layerPos += delta;
     
    136125      pred = pred - layerDist; // finally preds index
    137126
    138       while(data[pred].Key.CompareTo(bubble) > 0) {  // must terminate since inf at root
     127      while (data[pred].Key.CompareTo(bubble) > 0) {  // must terminate since inf at root
    139128        data[hole].Key = data[pred].Key;
    140129        data[hole].Value = data[pred].Value;
     
    149138      data[hole].Key = bubble;
    150139      data[hole].Value = data[sz].Value;
    151 
    152140      data[sz].Key = Supremum; // mark as deleted
    153141    }
     
    158146      finalLayerDist--;
    159147
    160       if(finalLayerDist == -1) { // layer full
    161                                  // start next layer
     148      if (finalLayerDist == -1) { // layer full
     149                                  // start next layer
    162150        finalLayerSize <<= 2;
    163151        finalLayerDist = finalLayerSize - 1;
     
    172160      pred = pred - layerDist; // finally preds index
    173161      var predKey = data[pred].Key;
    174       while(predKey.CompareTo(key) > 0) {
     162      while (predKey.CompareTo(key) > 0) {
    175163        data[hole].Key = predKey;
    176164        data[hole].Value = data[pred].Value;
     
    195183      var sup = Supremum;
    196184      var cap = capacity;
    197       for(var i = 1; i <= cap; i++) {
     185      for (var i = 1; i <= cap; i++) {
    198186        data[i].Key = sup;
    199187      }
  • trunk/sources/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/SpacePartitioningTree.cs

    r14862 r15207  
    5656using System;
    5757using System.Collections.Generic;
    58 using System.Linq;
    59 using HeuristicLab.Common;
    6058
    6159namespace HeuristicLab.Algorithms.DataAnalysis {
     
    6361  /// Space partitioning tree (SPTree)
    6462  /// </summary>
    65   public class SpacePartitioningTree : ISpacePartitioningTree {
    66     private const uint QT_NODE_CAPACITY = 1;
    67 
    68     private double[] buff;
    69     private SpacePartitioningTree parent;
     63  internal class SpacePartitioningTree {
     64    private const uint QtNodeCapacity = 1;
     65
     66    #region Fields
    7067    private int dimension;
    7168    private bool isLeaf;
     
    8077    // Indices in this space-partitioning tree node, corresponding center-of-mass, and list of all children
    8178    private double[] centerOfMass;
    82     private readonly int[] index = new int[QT_NODE_CAPACITY];
     79    private readonly int[] index = new int[QtNodeCapacity];
    8380
    8481    // Children
    8582    private SpacePartitioningTree[] children;
    8683    private uint noChildren;
     84    #endregion
    8785
    8886    public SpacePartitioningTree(double[,] inpData) {
     
    104102      var width = new double[d];
    105103      for (var i = 0; i < d; i++) width[i] = Math.Max(maxY[i] - meanY[i], meanY[i] - minY[i]) + 1e-5;
    106       Init(null, inpData, meanY, width);
     104      Init(inpData, meanY, width);
    107105      Fill(n);
    108106    }
    109107
    110     public SpacePartitioningTree(double[,] inpData, IEnumerable<double> impCorner, IEnumerable<double> impWith) {
    111       Init(null, inpData, impCorner, impWith);
    112     }
    113     public SpacePartitioningTree(SpacePartitioningTree parent, double[,] inpData, IEnumerable<double> impCorner, IEnumerable<double> impWith) {
    114       Init(parent, inpData, impCorner, impWith);
    115     }
    116 
    117     public ISpacePartitioningTree GetParent() {
    118       return parent;
     108    private SpacePartitioningTree(double[,] inpData, IEnumerable<double> impCorner, IEnumerable<double> impWith) {
     109      Init(inpData, impCorner, impWith);
    119110    }
    120111
     
    132123
    133124      // If there is space in this quad tree and it is a leaf, add the object here
    134       if (isLeaf && size < QT_NODE_CAPACITY) {
     125      if (isLeaf && size < QtNodeCapacity) {
    135126        index[size] = newIndex;
    136127        size++;
     
    138129      }
    139130
    140       // Don't add duplicates for now (this is not very nice)
     131      // Don't add duplicates
    141132      var anyDuplicate = false;
    142133      for (uint n = 0; n < size; n++) {
     
    161152    }
    162153
    163     public void Subdivide() {
    164       // Create new children
    165       var newCorner = new double[dimension];
    166       var newWidth = new double[dimension];
    167       for (var i = 0; i < noChildren; i++) {
    168         var div = 1;
    169         for (var d = 0; d < dimension; d++) {
    170           newWidth[d] = .5 * boundary.GetWidth(d);
    171           if ((i / div) % 2 == 1) newCorner[d] = boundary.GetCorner(d) - .5 * boundary.GetWidth(d);
    172           else newCorner[d] = boundary.GetCorner(d) + .5 * boundary.GetWidth(d);
    173           div *= 2;
    174         }
    175         children[i] = new SpacePartitioningTree(this, data, newCorner, newWidth);
    176       }
    177 
    178       // Move existing points to correct children
    179       for (var i = 0; i < size; i++) {
    180         var success = false;
    181         for (var j = 0; j < noChildren; j++) {
    182           if (!success) success = children[j].Insert(index[i]);
    183         }
    184         index[i] = -1; // as in tSNE implementation by van der Maaten
    185       }
    186 
    187       // Empty parent node
    188       size = 0;
    189       isLeaf = false;
    190     }
    191 
    192     public bool IsCorrect() {
    193       var row = new double[dimension];
    194       for (var n = 0; n < size; n++) {
    195         Buffer.BlockCopy(data, sizeof(double) * dimension * index[n], row, 0, sizeof(double) * dimension);
    196         if (!boundary.ContainsPoint(row)) return false;
    197       }
    198       if (isLeaf) return true;
    199       var correct = true;
    200       for (var i = 0; i < noChildren; i++) correct = correct && children[i].IsCorrect();
    201       return correct;
    202     }
    203 
    204     public void GetAllIndices(int[] indices) {
    205       GetAllIndices(indices, 0);
    206     }
    207 
    208     public int GetAllIndices(int[] indices, int loc) {
    209       // Gather indices in current quadrant
    210       for (var i = 0; i < size; i++) indices[loc + i] = index[i];
    211       loc += (int)size;
    212       // Gather indices in children
    213       if (isLeaf) return loc;
    214       for (var i = 0; i < noChildren; i++) loc = children[i].GetAllIndices(indices, loc);
    215       return loc;
    216     }
    217 
    218     public int GetDepth() {
    219       return isLeaf ? 1 : 1 + children.Max(x => x.GetDepth());
    220     }
    221 
    222154    public void ComputeNonEdgeForces(int pointIndex, double theta, double[] negF, ref double sumQ) {
    223155      // Make sure that we spend no time on empty nodes or self-interactions
     
    226158      // Compute distance between point and center-of-mass
    227159      var D = .0;
     160      var buff = new double[dimension];
    228161      for (var d = 0; d < dimension; d++) buff[d] = data[pointIndex, d] - centerOfMass[d];
    229162      for (var d = 0; d < dimension; d++) D += buff[d] * buff[d];
     
    233166      for (var d = 0; d < dimension; d++) {
    234167        var curWidth = boundary.GetWidth(d);
    235         maxWidth = (maxWidth > curWidth) ? maxWidth : curWidth;
     168        maxWidth = maxWidth > curWidth ? maxWidth : curWidth;
    236169      }
    237170      if (isLeaf || maxWidth / Math.Sqrt(D) < theta) {
     
    250183    }
    251184
    252     // does not use the tree
    253     public void ComputeEdgeForces(int[] rowP, int[] colP, double[] valP, int n, double[,] posF) {
     185    public static void ComputeEdgeForces(int[] rowP, int[] colP, double[] valP, int n, double[,] posF, double[,] data, int dimension) {
    254186      // Loop over all edges in the graph
    255187      for (var k = 0; k < n; k++) {
     
    259191          // uses squared distance
    260192          var d = 1.0;
     193          var buff = new double[dimension];
    261194          for (var j = 0; j < dimension; j++) buff[j] = data[k, j] - data[colP[i], j];
    262195          for (var j = 0; j < dimension; j++) d += buff[j] * buff[j];
     
    273206      for (var i = 0; i < n; i++) Insert(i);
    274207    }
    275     private void Init(SpacePartitioningTree p, double[,] inpData, IEnumerable<double> inpCorner, IEnumerable<double> inpWidth) {
    276       parent = p;
     208
     209    private void Init(double[,] inpData, IEnumerable<double> inpCorner, IEnumerable<double> inpWidth) {
    277210      dimension = inpData.GetLength(1);
    278211      noChildren = 2;
     
    283216      cumulativeSize = 0;
    284217      boundary = new Cell((uint)dimension);
     218
    285219      inpCorner.ForEach((i, x) => boundary.SetCorner(i, x));
    286220      inpWidth.ForEach((i, x) => boundary.SetWidth(i, x));
     
    288222      children = new SpacePartitioningTree[noChildren];
    289223      centerOfMass = new double[dimension];
    290       buff = new double[dimension];
    291 
     224    }
     225
     226    private void Subdivide() {
     227      // Create new children
     228      var newCorner = new double[dimension];
     229      var newWidth = new double[dimension];
     230      for (var i = 0; i < noChildren; i++) {
     231        var div = 1;
     232        for (var d = 0; d < dimension; d++) {
     233          newWidth[d] = .5 * boundary.GetWidth(d);
     234          if (i / div % 2 == 1) newCorner[d] = boundary.GetCorner(d) - .5 * boundary.GetWidth(d);
     235          else newCorner[d] = boundary.GetCorner(d) + .5 * boundary.GetWidth(d);
     236          div *= 2;
     237        }
     238        children[i] = new SpacePartitioningTree(data, newCorner, newWidth);
     239      }
     240
     241      // Move existing points to correct children
     242      for (var i = 0; i < size; i++) {
     243        var success = false;
     244        for (var j = 0; j < noChildren; j++) {
     245          if (!success) success = children[j].Insert(index[i]);
     246        }
     247        index[i] = -1; // as in tSNE implementation by van der Maaten
     248      }
     249      // Empty parent node
     250      size = 0;
     251      isLeaf = false;
    292252    }
    293253    #endregion
    294 
    295254
    296255    private class Cell {
  • trunk/sources/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/TSNEAlgorithm.cs

    r14863 r15207  
    4141  /// </summary>
    4242  [Item("tSNE", "t-distributed stochastic neighbourhood embedding projects the data in a low " +
    43                 "dimensional space to allow visual cluster identification.")]
     43                "dimensional space to allow visual cluster identification. Implemented similar to: https://lvdmaaten.github.io/tsne/#implementations (Barnes-Hut t-SNE). Described in : https://lvdmaaten.github.io/publications/papers/JMLR_2014.pdf")]
    4444  [Creatable(CreatableAttribute.Categories.DataAnalysis, Priority = 100)]
    4545  [StorableClass]
     
    198198
    199199    private TSNEAlgorithm(TSNEAlgorithm original, Cloner cloner) : base(original, cloner) {
    200       if(original.dataRowNames!=null)
    201       this.dataRowNames = new Dictionary<string, List<int>>(original.dataRowNames);
     200      if (original.dataRowNames != null)
     201        this.dataRowNames = new Dictionary<string, List<int>>(original.dataRowNames);
    202202      if (original.dataRows != null)
    203203        this.dataRows = original.dataRows.ToDictionary(kvp => kvp.Key, kvp => cloner.Clone(kvp.Value));
     
    226226      Parameters.Add(new FixedValueParameter<BoolValue>(SetSeedRandomlyParameterName, "If the seed should be random.", new BoolValue(true)));
    227227      Parameters.Add(new FixedValueParameter<IntValue>(SeedParameterName, "The seed used if it should not be random.", new IntValue(0)));
    228       Parameters.Add(new FixedValueParameter<StringValue>(ClassesParameterName, "name of the column specifying the class lables of each data point. If the label column can not be found training/test is used as labels.", new StringValue("none")));
     228      Parameters.Add(new FixedValueParameter<StringValue>(ClassesParameterName, "Name of the column specifying the class lables of each data point. If the label column can not be found training/test is used as labels.", new StringValue("none")));
    229229      Parameters.Add(new FixedValueParameter<BoolValue>(NormalizationParameterName, "Whether the data should be zero centered and have variance of 1 for each variable, so different scalings are ignored.", new BoolValue(true)));
    230230      Parameters.Add(new FixedValueParameter<IntValue>(UpdateIntervalParameterName, "", new IntValue(50)));
  • trunk/sources/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/TSNEStatic.cs

    r14862 r15207  
    5656using System;
    5757using System.Collections.Generic;
    58 using System.Linq;
    5958using HeuristicLab.Collections;
    6059using HeuristicLab.Common;
    6160using HeuristicLab.Core;
    62 using HeuristicLab.Optimization;
    6361using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
    6462using HeuristicLab.Random;
     
    7068    [StorableClass]
    7169    public sealed class TSNEState : DeepCloneable {
     70      #region Storables
    7271      // initialized once
    7372      [Storable]
     
    122121      [Storable]
    123122      public double[,] dY;
    124 
     123      #endregion
     124
     125      #region Constructors & Cloning
    125126      private TSNEState(TSNEState original, Cloner cloner) : base(original, cloner) {
    126         this.distance = cloner.Clone(original.distance);
    127         this.random = cloner.Clone(original.random);
    128         this.perplexity = original.perplexity;
    129         this.exact = original.exact;
    130         this.noDatapoints = original.noDatapoints;
    131         this.finalMomentum = original.finalMomentum;
    132         this.momSwitchIter = original.momSwitchIter;
    133         this.stopLyingIter = original.stopLyingIter;
    134         this.theta = original.theta;
    135         this.eta = original.eta;
    136         this.newDimensions = original.newDimensions;
    137         if(original.valP != null) {
    138           this.valP = new double[original.valP.Length];
    139           Array.Copy(original.valP, this.valP, this.valP.Length);
    140         }
    141         if(original.rowP != null) {
    142           this.rowP = new int[original.rowP.Length];
    143           Array.Copy(original.rowP, this.rowP, this.rowP.Length);
    144         }
    145         if(original.colP != null) {
    146           this.colP = new int[original.colP.Length];
    147           Array.Copy(original.colP, this.colP, this.colP.Length);
    148         }
    149         if(original.p != null) {
    150           this.p = new double[original.p.GetLength(0), original.p.GetLength(1)];
    151           Array.Copy(original.p, this.p, this.p.Length);
    152         }
    153         this.newData = new double[original.newData.GetLength(0), original.newData.GetLength(1)];
    154         Array.Copy(original.newData, this.newData, this.newData.Length);
    155         this.iter = original.iter;
    156         this.currentMomentum = original.currentMomentum;
    157         this.gains = new double[original.gains.GetLength(0), original.gains.GetLength(1)];
    158         Array.Copy(original.gains, this.gains, this.gains.Length);
    159         this.uY = new double[original.uY.GetLength(0), original.uY.GetLength(1)];
    160         Array.Copy(original.uY, this.uY, this.uY.Length);
    161         this.dY = new double[original.dY.GetLength(0), original.dY.GetLength(1)];
    162         Array.Copy(original.dY, this.dY, this.dY.Length);
     127        distance = cloner.Clone(original.distance);
     128        random = cloner.Clone(original.random);
     129        perplexity = original.perplexity;
     130        exact = original.exact;
     131        noDatapoints = original.noDatapoints;
     132        finalMomentum = original.finalMomentum;
     133        momSwitchIter = original.momSwitchIter;
     134        stopLyingIter = original.stopLyingIter;
     135        theta = original.theta;
     136        eta = original.eta;
     137        newDimensions = original.newDimensions;
     138        if (original.valP != null) {
     139          valP = new double[original.valP.Length];
     140          Array.Copy(original.valP, valP, valP.Length);
     141        }
     142        if (original.rowP != null) {
     143          rowP = new int[original.rowP.Length];
     144          Array.Copy(original.rowP, rowP, rowP.Length);
     145        }
     146        if (original.colP != null) {
     147          colP = new int[original.colP.Length];
     148          Array.Copy(original.colP, colP, colP.Length);
     149        }
     150        if (original.p != null) {
     151          p = new double[original.p.GetLength(0), original.p.GetLength(1)];
     152          Array.Copy(original.p, p, p.Length);
     153        }
     154        newData = new double[original.newData.GetLength(0), original.newData.GetLength(1)];
     155        Array.Copy(original.newData, newData, newData.Length);
     156        iter = original.iter;
     157        currentMomentum = original.currentMomentum;
     158        gains = new double[original.gains.GetLength(0), original.gains.GetLength(1)];
     159        Array.Copy(original.gains, gains, gains.Length);
     160        uY = new double[original.uY.GetLength(0), original.uY.GetLength(1)];
     161        Array.Copy(original.uY, uY, uY.Length);
     162        dY = new double[original.dY.GetLength(0), original.dY.GetLength(1)];
     163        Array.Copy(original.dY, dY, dY.Length);
    163164      }
    164165
     
    177178        this.stopLyingIter = stopLyingIter;
    178179        this.momSwitchIter = momSwitchIter;
    179         this.currentMomentum = momentum;
     180        currentMomentum = momentum;
    180181        this.finalMomentum = finalMomentum;
    181182        this.eta = eta;
    182183
    183 
    184184        // initialize
    185185        noDatapoints = data.Length;
    186         if(noDatapoints - 1 < 3 * perplexity)
     186        if (noDatapoints - 1 < 3 * perplexity)
    187187          throw new ArgumentException("Perplexity too large for the number of data points!");
    188188
     
    192192        uY = new double[noDatapoints, newDimensions];
    193193        gains = new double[noDatapoints, newDimensions];
    194         for(var i = 0; i < noDatapoints; i++)
    195           for(var j = 0; j < newDimensions; j++)
     194        for (var i = 0; i < noDatapoints; i++)
     195          for (var j = 0; j < newDimensions; j++)
    196196            gains[i, j] = 1.0;
    197197
     
    206206
    207207        // Lie about the P-values (factor is 4 in the MATLAB implementation)
    208         if(exact) for(var i = 0; i < noDatapoints; i++) for(var j = 0; j < noDatapoints; j++) p[i, j] *= 12.0;
    209         else for(var i = 0; i < rowP[noDatapoints]; i++) valP[i] *= 12.0;
     208        if (exact) for (var i = 0; i < noDatapoints; i++) for (var j = 0; j < noDatapoints; j++) p[i, j] *= 12.0;
     209        else for (var i = 0; i < rowP[noDatapoints]; i++) valP[i] *= 12.0;
    210210
    211211        // Initialize solution (randomly)
    212212        var rand = new NormalDistributedRandom(random, 0, 1);
    213         for(var i = 0; i < noDatapoints; i++)
    214           for(var j = 0; j < newDimensions; j++)
     213        for (var i = 0; i < noDatapoints; i++)
     214          for (var j = 0; j < newDimensions; j++)
    215215            newData[i, j] = rand.NextDouble() * .0001;
    216216      }
     217      #endregion
    217218
    218219      public double EvaluateError() {
     
    222223      }
    223224
     225      #region Helpers
    224226      private static void CalculateApproximateSimilarities(T[] data, IDistance<T> distance, double perplexity, out int[] rowP, out int[] colP, out double[] valP) {
    225227        // Compute asymmetric pairwise input similarities
     
    233235        valP = sValP;
    234236        var sumP = .0;
    235         for(var i = 0; i < rowP[data.Length]; i++) sumP += valP[i];
    236         for(var i = 0; i < rowP[data.Length]; i++) valP[i] /= sumP;
     237        for (var i = 0; i < rowP[data.Length]; i++) sumP += valP[i];
     238        for (var i = 0; i < rowP[data.Length]; i++) valP[i] /= sumP;
    237239      }
    238240
     
    242244        ComputeGaussianPerplexity(data, distance, p, perplexity);
    243245        // Symmetrize input similarities
    244         for(var n = 0; n < data.Length; n++) {
    245           for(var m = n + 1; m < data.Length; m++) {
     246        for (var n = 0; n < data.Length; n++) {
     247          for (var m = n + 1; m < data.Length; m++) {
    246248            p[n, m] += p[m, n];
    247249            p[m, n] = p[n, m];
     
    249251        }
    250252        var sumP = .0;
    251         for(var i = 0; i < data.Length; i++) for(var j = 0; j < data.Length; j++) sumP += p[i, j];
    252         for(var i = 0; i < data.Length; i++) for(var j = 0; j < data.Length; j++) p[i, j] /= sumP;
     253        for (var i = 0; i < data.Length; i++) for (var j = 0; j < data.Length; j++) sumP += p[i, j];
     254        for (var i = 0; i < data.Length; i++) for (var j = 0; j < data.Length; j++) p[i, j] /= sumP;
    253255        return p;
    254256      }
    255257
    256258      private static void ComputeGaussianPerplexity(IReadOnlyList<T> x, IDistance<T> distance, out int[] rowP, out int[] colP, out double[] valP, double perplexity, int k) {
    257         if(perplexity > k) throw new ArgumentException("Perplexity should be lower than k!");
    258 
    259         int n = x.Count;
     259        if (perplexity > k) throw new ArgumentException("Perplexity should be lower than k!");
     260
     261        var n = x.Count;
    260262        // Allocate the memory we need
    261263        rowP = new int[n + 1];
     
    264266        var curP = new double[n - 1];
    265267        rowP[0] = 0;
    266         for(var i = 0; i < n; i++) rowP[i + 1] = rowP[i] + k;
     268        for (var i = 0; i < n; i++) rowP[i + 1] = rowP[i] + k;
    267269
    268270        var objX = new List<IndexedItem<T>>();
    269         for(var i = 0; i < n; i++) objX.Add(new IndexedItem<T>(i, x[i]));
     271        for (var i = 0; i < n; i++) objX.Add(new IndexedItem<T>(i, x[i]));
    270272
    271273        // Build ball tree on data set
     
    273275
    274276        // Loop over all points to find nearest neighbors
    275         for(var i = 0; i < n; i++) {
     277        for (var i = 0; i < n; i++) {
    276278          IList<IndexedItem<T>> indices;
    277279          IList<double> distances;
     
    285287          var minBeta = double.MinValue;
    286288          var maxBeta = double.MaxValue;
    287           const double tol = 1e-5; 
     289          const double tol = 1e-5;
    288290
    289291          // Iterate until we found a good perplexity
    290292          var iter = 0; double sumP = 0;
    291           while(!found && iter < 200) { 
     293          while (!found && iter < 200) {
    292294
    293295            // Compute Gaussian kernel row
    294             for(var m = 0; m < k; m++) curP[m] = Math.Exp(-beta * distances[m + 1]);
     296            for (var m = 0; m < k; m++) curP[m] = Math.Exp(-beta * distances[m + 1]);
    295297
    296298            // Compute entropy of current row
    297299            sumP = double.Epsilon;
    298             for(var m = 0; m < k; m++) sumP += curP[m];
     300            for (var m = 0; m < k; m++) sumP += curP[m];
    299301            var h = .0;
    300             for(var m = 0; m < k; m++) h += beta * (distances[m + 1] * curP[m]);
     302            for (var m = 0; m < k; m++) h += beta * (distances[m + 1] * curP[m]);
    301303            h = h / sumP + Math.Log(sumP);
    302304
    303305            // Evaluate whether the entropy is within the tolerance level
    304306            var hdiff = h - Math.Log(perplexity);
    305             if(hdiff < tol && -hdiff < tol) {
     307            if (hdiff < tol && -hdiff < tol) {
    306308              found = true;
    307309            } else {
    308               if(hdiff > 0) {
     310              if (hdiff > 0) {
    309311                minBeta = beta;
    310                 if(maxBeta.IsAlmost(double.MaxValue) || maxBeta.IsAlmost(double.MinValue))
     312                if (maxBeta.IsAlmost(double.MaxValue) || maxBeta.IsAlmost(double.MinValue))
    311313                  beta *= 2.0;
    312314                else
     
    314316              } else {
    315317                maxBeta = beta;
    316                 if(minBeta.IsAlmost(double.MinValue) || minBeta.IsAlmost(double.MaxValue))
     318                if (minBeta.IsAlmost(double.MinValue) || minBeta.IsAlmost(double.MaxValue))
    317319                  beta /= 2.0;
    318320                else
     
    326328
    327329          // Row-normalize current row of P and store in matrix
    328           for(var m = 0; m < k; m++) curP[m] /= sumP;
    329           for(var m = 0; m < k; m++) {
     330          for (var m = 0; m < k; m++) curP[m] /= sumP;
     331          for (var m = 0; m < k; m++) {
    330332            colP[rowP[i] + m] = indices[m + 1].Index;
    331333            valP[rowP[i] + m] = curP[m];
     
    337339        var dd = ComputeDistances(x, distance);
    338340
    339         int n = x.Length;
     341        var n = x.Length;
    340342        // Compute the Gaussian kernel row by row
    341         for(var i = 0; i < n; i++) {
     343        for (var i = 0; i < n; i++) {
    342344          // Initialize some variables
    343345          var found = false;
     
    350352          // Iterate until we found a good perplexity
    351353          var iter = 0;
    352           while(!found && iter < 200) {      // 200 iterations as in tSNE implementation by van der Maarten
     354          while (!found && iter < 200) {      // 200 iterations as in tSNE implementation by van der Maarten
    353355
    354356            // Compute Gaussian kernel row
    355             for(var m = 0; m < n; m++) p[i, m] = Math.Exp(-beta * dd[i][m]);
     357            for (var m = 0; m < n; m++) p[i, m] = Math.Exp(-beta * dd[i][m]);
    356358            p[i, i] = double.Epsilon;
    357359
    358360            // Compute entropy of current row
    359361            sumP = double.Epsilon;
    360             for(var m = 0; m < n; m++) sumP += p[i, m];
     362            for (var m = 0; m < n; m++) sumP += p[i, m];
    361363            var h = 0.0;
    362             for(var m = 0; m < n; m++) h += beta * (dd[i][m] * p[i, m]);
     364            for (var m = 0; m < n; m++) h += beta * (dd[i][m] * p[i, m]);
    363365            h = h / sumP + Math.Log(sumP);
    364366
    365367            // Evaluate whether the entropy is within the tolerance level
    366368            var hdiff = h - Math.Log(perplexity);
    367             if(hdiff < tol && -hdiff < tol) {
     369            if (hdiff < tol && -hdiff < tol) {
    368370              found = true;
    369371            } else {
    370               if(hdiff > 0) {
     372              if (hdiff > 0) {
    371373                minBeta = beta;
    372                 if(maxBeta.IsAlmost(double.MaxValue) || maxBeta.IsAlmost(double.MinValue))
     374                if (maxBeta.IsAlmost(double.MaxValue) || maxBeta.IsAlmost(double.MinValue))
    373375                  beta *= 2.0;
    374376                else
     
    376378              } else {
    377379                maxBeta = beta;
    378                 if(minBeta.IsAlmost(double.MinValue) || minBeta.IsAlmost(double.MaxValue))
     380                if (minBeta.IsAlmost(double.MinValue) || minBeta.IsAlmost(double.MaxValue))
    379381                  beta /= 2.0;
    380382                else
     
    388390
    389391          // Row normalize P
    390           for(var m = 0; m < n; m++) p[i, m] /= sumP;
     392          for (var m = 0; m < n; m++) p[i, m] /= sumP;
    391393        }
    392394      }
     
    394396      private static double[][] ComputeDistances(T[] x, IDistance<T> distance) {
    395397        var res = new double[x.Length][];
    396         for(int r = 0; r < x.Length; r++) {
     398        for (var r = 0; r < x.Length; r++) {
    397399          var rowV = new double[x.Length];
    398400          // all distances must be symmetric
    399           for(int c = 0; c < r; c++) {
     401          for (var c = 0; c < r; c++) {
    400402            rowV[c] = res[c][r];
    401403          }
    402404          rowV[r] = 0.0; // distance to self is zero for all distances
    403           for(int c = r + 1; c < x.Length; c++) {
     405          for (var c = r + 1; c < x.Length; c++) {
    404406            rowV[c] = distance.Get(x[r], x[c]);
    405407          }
     
    418420        // Compute Q-matrix and normalization sum
    419421        var sumQ = double.Epsilon;
    420         for(var n1 = 0; n1 < n; n1++) {
    421           for(var m = 0; m < n; m++) {
    422             if(n1 != m) {
     422        for (var n1 = 0; n1 < n; n1++) {
     423          for (var m = 0; m < n; m++) {
     424            if (n1 != m) {
    423425              q[n1, m] = 1 / (1 + dd[n1, m]);
    424426              sumQ += q[n1, m];
     
    426428          }
    427429        }
    428         for(var i = 0; i < n; i++) for(var j = 0; j < n; j++) q[i, j] /= sumQ;
     430        for (var i = 0; i < n; i++) for (var j = 0; j < n; j++) q[i, j] /= sumQ;
    429431
    430432        // Sum t-SNE error
    431433        var c = .0;
    432         for(var i = 0; i < n; i++)
    433           for(var j = 0; j < n; j++) {
     434        for (var i = 0; i < n; i++)
     435          for (var j = 0; j < n; j++) {
    434436            c += p[i, j] * Math.Log((p[i, j] + float.Epsilon) / (q[i, j] + float.Epsilon));
    435437          }
     
    437439      }
    438440
    439       // The mapping of the approximate tSNE looks good but the error curve never changes.
    440441      private static double EvaluateErrorApproximate(IReadOnlyList<int> rowP, IReadOnlyList<int> colP, IReadOnlyList<double> valP, double[,] y, double theta) {
    441442        // Get estimate of normalization term
     
    444445        var tree = new SpacePartitioningTree(y);
    445446        var buff = new double[d];
    446         double sumQ = 0.0;
    447         for(var i = 0; i < n; i++) tree.ComputeNonEdgeForces(i, theta, buff, ref sumQ);
     447        var sumQ = 0.0;
     448        for (var i = 0; i < n; i++) tree.ComputeNonEdgeForces(i, theta, buff, ref sumQ);
    448449
    449450        // Loop over all edges to compute t-SNE error
    450451        var c = .0;
    451         for(var k = 0; k < n; k++) {
    452           for(var i = rowP[k]; i < rowP[k + 1]; i++) {
     452        for (var k = 0; k < n; k++) {
     453          for (var i = rowP[k]; i < rowP[k + 1]; i++) {
    453454            var q = .0;
    454             for(var j = 0; j < d; j++) buff[j] = y[k, j];
    455             for(var j = 0; j < d; j++) buff[j] -= y[colP[i], j];
    456             for(var j = 0; j < d; j++) q += buff[j] * buff[j];
     455            for (var j = 0; j < d; j++) buff[j] = y[k, j];
     456            for (var j = 0; j < d; j++) buff[j] -= y[colP[i], j];
     457            for (var j = 0; j < d; j++) q += buff[j] * buff[j];
    457458            q = (1.0 / (1.0 + q)) / sumQ;
    458459            c += valP[i] * Math.Log((valP[i] + float.Epsilon) / (q + float.Epsilon));
     
    466467        var n = rowP.Count - 1;
    467468        var rowCounts = new int[n];
    468         for(var j = 0; j < n; j++) {
    469           for(var i = rowP[j]; i < rowP[j + 1]; i++) {
     469        for (var j = 0; j < n; j++) {
     470          for (var i = rowP[j]; i < rowP[j + 1]; i++) {
    470471
    471472            // Check whether element (col_P[i], n) is present
    472473            var present = false;
    473             for(var m = rowP[colP[i]]; m < rowP[colP[i] + 1]; m++) {
    474               if(colP[m] == j) present = true;
     474            for (var m = rowP[colP[i]]; m < rowP[colP[i] + 1]; m++) {
     475              if (colP[m] == j) present = true;
    475476            }
    476             if(present) rowCounts[j]++;
     477            if (present) rowCounts[j]++;
    477478            else {
    478479              rowCounts[j]++;
     
    482483        }
    483484        var noElem = 0;
    484         for(var i = 0; i < n; i++) noElem += rowCounts[i];
     485        for (var i = 0; i < n; i++) noElem += rowCounts[i];
    485486
    486487        // Allocate memory for symmetrized matrix
     
    491492        // Construct new row indices for symmetric matrix
    492493        symRowP[0] = 0;
    493         for(var i = 0; i < n; i++) symRowP[i + 1] = symRowP[i] + rowCounts[i];
     494        for (var i = 0; i < n; i++) symRowP[i + 1] = symRowP[i] + rowCounts[i];
    494495
    495496        // Fill the result matrix
    496497        var offset = new int[n];
    497         for(var j = 0; j < n; j++) {
    498           for(var i = rowP[j]; i < rowP[j + 1]; i++) {                                  // considering element(n, colP[i])
     498        for (var j = 0; j < n; j++) {
     499          for (var i = rowP[j]; i < rowP[j + 1]; i++) {                                  // considering element(n, colP[i])
    499500
    500501            // Check whether element (col_P[i], n) is present
    501502            var present = false;
    502             for(var m = rowP[colP[i]]; m < rowP[colP[i] + 1]; m++) {
    503               if(colP[m] != j) continue;
     503            for (var m = rowP[colP[i]]; m < rowP[colP[i] + 1]; m++) {
     504              if (colP[m] != j) continue;
    504505              present = true;
    505               if(j > colP[i]) continue; // make sure we do not add elements twice
     506              if (j > colP[i]) continue; // make sure we do not add elements twice
    506507              symColP[symRowP[j] + offset[j]] = colP[i];
    507508              symColP[symRowP[colP[i]] + offset[colP[i]]] = j;
     
    511512
    512513            // If (colP[i], n) is not present, there is no addition involved
    513             if(!present) {
     514            if (!present) {
    514515              symColP[symRowP[j] + offset[j]] = colP[i];
    515516              symColP[symRowP[colP[i]] + offset[colP[i]]] = j;
     
    519520
    520521            // Update offsets
    521             if(present && (j > colP[i])) continue;
     522            if (present && (j > colP[i])) continue;
    522523            offset[j]++;
    523             if(colP[i] != j) offset[colP[i]]++;
    524           }
    525         }
    526 
    527         for(var i = 0; i < noElem; i++) symValP[i] /= 2.0;
    528       }
     524            if (colP[i] != j) offset[colP[i]]++;
     525          }
     526        }
     527
     528        for (var i = 0; i < noElem; i++) symValP[i] /= 2.0;
     529      }
     530      #endregion
    529531    }
    530532
    531533    /// <summary>
    532     /// Simple interface to tSNE
     534    /// Static interface to tSNE
    533535    /// </summary>
    534536    /// <param name="data"></param>
     
    548550      int newDimensions = 2, double perplexity = 25, int iterations = 1000,
    549551      double theta = 0,
    550       int stopLyingIter = 250, int momSwitchIter = 250, double momentum = .5,
    551       double finalMomentum = .8, double eta = 200.0
     552      int stopLyingIter = 0, int momSwitchIter = 0, double momentum = .5,
     553      double finalMomentum = .8, double eta = 10.0
    552554      ) {
    553555      var state = CreateState(data, distance, random, newDimensions, perplexity,
    554556        theta, stopLyingIter, momSwitchIter, momentum, finalMomentum, eta);
    555557
    556       for(int i = 0; i < iterations - 1; i++) {
     558      for (var i = 0; i < iterations - 1; i++) {
    557559        Iterate(state);
    558560      }
     
    562564    public static TSNEState CreateState(T[] data, IDistance<T> distance, IRandom random,
    563565      int newDimensions = 2, double perplexity = 25, double theta = 0,
    564       int stopLyingIter = 250, int momSwitchIter = 250, double momentum = .5,
    565       double finalMomentum = .8, double eta = 200.0
     566      int stopLyingIter = 0, int momSwitchIter = 0, double momentum = .5,
     567      double finalMomentum = .8, double eta = 10.0
    566568      ) {
    567569      return new TSNEState(data, distance, random, newDimensions, perplexity, theta, stopLyingIter, momSwitchIter, momentum, finalMomentum, eta);
    568570    }
    569571
    570 
    571572    public static double[,] Iterate(TSNEState state) {
    572       if(state.exact)
     573      if (state.exact)
    573574        ComputeExactGradient(state.p, state.newData, state.noDatapoints, state.newDimensions, state.dY);
    574575      else
     
    576577
    577578      // Update gains
    578       for(var i = 0; i < state.noDatapoints; i++) {
    579         for(var j = 0; j < state.newDimensions; j++) {
     579      for (var i = 0; i < state.noDatapoints; i++) {
     580        for (var j = 0; j < state.newDimensions; j++) {
    580581          state.gains[i, j] = Math.Sign(state.dY[i, j]) != Math.Sign(state.uY[i, j])
    581582            ? state.gains[i, j] + .2  // +0.2 nd *0.8 are used in two separate implementations of tSNE -> seems to be correct
    582583            : state.gains[i, j] * .8;
    583584
    584           if(state.gains[i, j] < .01) state.gains[i, j] = .01;
     585          if (state.gains[i, j] < .01) state.gains[i, j] = .01;
    585586        }
    586587      }
     
    588589
    589590      // Perform gradient update (with momentum and gains)
    590       for(var i = 0; i < state.noDatapoints; i++)
    591         for(var j = 0; j < state.newDimensions; j++)
     591      for (var i = 0; i < state.noDatapoints; i++)
     592        for (var j = 0; j < state.newDimensions; j++)
    592593          state.uY[i, j] = state.currentMomentum * state.uY[i, j] - state.eta * state.gains[i, j] * state.dY[i, j];
    593594
    594       for(var i = 0; i < state.noDatapoints; i++)
    595         for(var j = 0; j < state.newDimensions; j++)
     595      for (var i = 0; i < state.noDatapoints; i++)
     596        for (var j = 0; j < state.newDimensions; j++)
    596597          state.newData[i, j] = state.newData[i, j] + state.uY[i, j];
    597598
     
    600601
    601602      // Stop lying about the P-values after a while, and switch momentum
    602       if(state.iter == state.stopLyingIter) {
    603         if(state.exact)
    604           for(var i = 0; i < state.noDatapoints; i++)
    605             for(var j = 0; j < state.noDatapoints; j++)
     603      if (state.iter == state.stopLyingIter) {
     604        if (state.exact)
     605          for (var i = 0; i < state.noDatapoints; i++)
     606            for (var j = 0; j < state.noDatapoints; j++)
    606607              state.p[i, j] /= 12.0;
    607608        else
    608           for(var i = 0; i < state.rowP[state.noDatapoints]; i++)
     609          for (var i = 0; i < state.rowP[state.noDatapoints]; i++)
    609610            state.valP[i] /= 12.0;
    610611      }
    611612
    612       if(state.iter == state.momSwitchIter)
     613      if (state.iter == state.momSwitchIter)
    613614        state.currentMomentum = state.finalMomentum;
    614615
     
    617618    }
    618619
    619 
    620 
     620    #region Helpers
    621621    private static void ComputeApproximateGradient(int[] rowP, int[] colP, double[] valP, double[,] y, int n, int d, double[,] dC, double theta) {
    622622      var tree = new SpacePartitioningTree(y);
    623       double sumQ = 0.0;
     623      var sumQ = 0.0;
    624624      var posF = new double[n, d];
    625625      var negF = new double[n, d];
    626       tree.ComputeEdgeForces(rowP, colP, valP, n, posF);
     626      SpacePartitioningTree.ComputeEdgeForces(rowP, colP, valP, n, posF, y, d);
    627627      var row = new double[d];
    628       for(var n1 = 0; n1 < n; n1++) {
    629         Array.Clear(row,0,row.Length);
     628      for (var n1 = 0; n1 < n; n1++) {
     629        Array.Clear(row, 0, row.Length);
    630630        tree.ComputeNonEdgeForces(n1, theta, row, ref sumQ);
    631         Buffer.BlockCopy(row, 0, negF, (sizeof(double) * n1 * d), d*sizeof(double));
     631        Buffer.BlockCopy(row, 0, negF, (sizeof(double) * n1 * d), d * sizeof(double));
    632632      }
    633633
    634634      // Compute final t-SNE gradient
    635635      for (var i = 0; i < n; i++)
    636         for(var j = 0; j < d; j++) {
     636        for (var j = 0; j < d; j++) {
    637637          dC[i, j] = posF[i, j] - negF[i, j] / sumQ;
    638638        }
     
    640640
    641641    private static void ComputeExactGradient(double[,] p, double[,] y, int n, int d, double[,] dC) {
    642 
    643642      // Make sure the current gradient contains zeros
    644       for(var i = 0; i < n; i++) for(var j = 0; j < d; j++) dC[i, j] = 0.0;
     643      for (var i = 0; i < n; i++) for (var j = 0; j < d; j++) dC[i, j] = 0.0;
    645644
    646645      // Compute the squared Euclidean distance matrix
     
    651650      var q = new double[n, n];
    652651      var sumQ = .0;
    653       for(var n1 = 0; n1 < n; n1++) {
    654         for(var m = 0; m < n; m++) {
    655           if(n1 == m) continue;
     652      for (var n1 = 0; n1 < n; n1++) {
     653        for (var m = 0; m < n; m++) {
     654          if (n1 == m) continue;
    656655          q[n1, m] = 1 / (1 + dd[n1, m]);
    657656          sumQ += q[n1, m];
     
    660659
    661660      // Perform the computation of the gradient
    662       for(var n1 = 0; n1 < n; n1++) {
    663         for(var m = 0; m < n; m++) {
    664           if(n1 == m) continue;
     661      for (var n1 = 0; n1 < n; n1++) {
     662        for (var m = 0; m < n; m++) {
     663          if (n1 == m) continue;
    665664          var mult = (p[n1, m] - q[n1, m] / sumQ) * q[n1, m];
    666           for(var d1 = 0; d1 < d; d1++) {
     665          for (var d1 = 0; d1 < d; d1++) {
    667666            dC[n1, d1] += (y[n1, d1] - y[m, d1]) * mult;
    668667          }
     
    673672    private static void ComputeSquaredEuclideanDistance(double[,] x, int n, int d, double[,] dd) {
    674673      var dataSums = new double[n];
    675       for(var i = 0; i < n; i++) {
    676         for(var j = 0; j < d; j++) {
     674      for (var i = 0; i < n; i++) {
     675        for (var j = 0; j < d; j++) {
    677676          dataSums[i] += x[i, j] * x[i, j];
    678677        }
    679678      }
    680       for(var i = 0; i < n; i++) {
    681         for(var m = 0; m < n; m++) {
     679      for (var i = 0; i < n; i++) {
     680        for (var m = 0; m < n; m++) {
    682681          dd[i, m] = dataSums[i] + dataSums[m];
    683682        }
    684683      }
    685       for(var i = 0; i < n; i++) {
     684      for (var i = 0; i < n; i++) {
    686685        dd[i, i] = 0.0;
    687         for(var m = i + 1; m < n; m++) {
     686        for (var m = i + 1; m < n; m++) {
    688687          dd[i, m] = 0.0;
    689           for(var j = 0; j < d; j++) {
     688          for (var j = 0; j < d; j++) {
    690689            dd[i, m] += (x[i, j] - x[m, j]) * (x[i, j] - x[m, j]);
    691690          }
     
    700699      var d = x.GetLength(1);
    701700      var mean = new double[d];
    702       for(var i = 0; i < n; i++) {
    703         for(var j = 0; j < d; j++) {
     701      for (var i = 0; i < n; i++) {
     702        for (var j = 0; j < d; j++) {
    704703          mean[j] += x[i, j];
    705704        }
    706705      }
    707       for(var i = 0; i < d; i++) {
     706      for (var i = 0; i < d; i++) {
    708707        mean[i] /= n;
    709708      }
    710709      // Subtract data mean
    711       for(var i = 0; i < n; i++) {
    712         for(var j = 0; j < d; j++) {
     710      for (var i = 0; i < n; i++) {
     711        for (var j = 0; j < d; j++) {
    713712          x[i, j] -= mean[j];
    714713        }
    715714      }
    716715    }
     716    #endregion
    717717  }
    718718}
  • trunk/sources/HeuristicLab.Algorithms.DataAnalysis/3.4/TSNE/VantagePointTree.cs

    r14862 r15207  
    6868  /// </summary>
    6969  /// <typeparam name="T"></typeparam>
    70   public class VantagePointTree<T> : IVantagePointTree<T> {
    71     #region properties
     70  public class VantagePointTree<T> {
    7271    private readonly List<T> items;
    7372    private readonly Node root;
    7473    private readonly IDistance<T> distance;
    75     #endregion
    7674
    7775    public VantagePointTree(IDistance<T> distance, IEnumerable<T> items) : base() {
     
    8280    }
    8381
     82    /// <summary>
     83    /// provides the k-nearest neighbours to a certain target element
     84    /// </summary>
     85    /// <param name="target">The target element</param>
     86    /// <param name="k">How many neighbours</param>
     87    /// <param name="results">The nearest neighbouring elements</param>
     88    /// <param name="distances">The distances form the target corresponding to the neighbouring elements</param>
    8489    public void Search(T target, int k, out IList<T> results, out IList<double> distances) {
    8590      var heap = new PriorityQueue<double, IndexedItem<double>>(double.MaxValue, double.MinValue, k);
    86       double tau = double.MaxValue;
     91      var tau = double.MaxValue;
    8792      Search(root, target, k, heap, ref tau);
    8893      var res = new List<T>();
    8994      var dist = new List<double>();
    9095      while (heap.Size > 0) {
    91         res.Add(items[heap.PeekMinValue().Index]);          // actually max distance
     96        res.Add(items[heap.PeekMinValue().Index]);// actually max distance
    9297        dist.Add(heap.PeekMinValue().Value);
    9398        heap.RemoveMin();
    9499      }
    95       res.Reverse(); 
     100      res.Reverse();
    96101      dist.Reverse();
    97102      results = res;
     
    104109      if (dist < tau) {
    105110        if (heap.Size == k) heap.RemoveMin();   // remove furthest node from result list (if we already have k results)
    106         heap.Insert(-dist, new IndexedItem<double>(node.index, dist));     
     111        heap.Insert(-dist, new IndexedItem<double>(node.index, dist));
    107112        if (heap.Size == k) tau = heap.PeekMinValue().Value;
    108113      }
Note: See TracChangeset for help on using the changeset viewer.