From a56b0618bf7f91fddbcdb4ed702fb4dcc31a139d Mon Sep 17 00:00:00 2001
From: Sergey Lisitsyn <lisitsyn.s.o@gmail.com>
Date: Fri, 24 May 2013 00:45:39 +0400
Subject: [PATCH] Fixes for tapkee for proper external random shuffle and
 randomized eigendecomposition

---
 src/shogun/lib/tapkee/defines/random.hpp      | 25 ++++++++++---------
 src/shogun/lib/tapkee/methods.hpp             |  6 ++---
 .../tapkee/routines/eigendecomposition.hpp    |  3 +--
 .../tapkee/routines/manifold_sculpting.hpp    | 15 +++++------
 src/shogun/lib/tapkee/utils/features.hpp      | 11 ++++----
 5 files changed, 31 insertions(+), 29 deletions(-)

diff --git a/src/shogun/lib/tapkee/defines/random.hpp b/src/shogun/lib/tapkee/defines/random.hpp
index 3419b3e716b..59ffa10b1f9 100644
--- a/src/shogun/lib/tapkee/defines/random.hpp
+++ b/src/shogun/lib/tapkee/defines/random.hpp
@@ -8,45 +8,46 @@
 
 #include <cstdlib>
 #include <algorithm>
+#include <limits>
 
 namespace tapkee
 {
 
-IndexType uniform_random_index() 
+inline IndexType uniform_random_index() 
 {
 #ifdef CUSTOM_UNIFORM_RANDOM_INDEX_FUNCTION
-	return CUSTOM_UNIFORM_RANDOM_INDEX_FUNCTION;
+	return CUSTOM_UNIFORM_RANDOM_INDEX_FUNCTION % std::numeric_limits<IndexType>::max();
 #else
-	return rand();
+	return std::rand();
 #endif
 }
 
-IndexType uniform_random_index_bounded(IndexType upper)
+inline IndexType uniform_random_index_bounded(IndexType upper)
 {
 	return uniform_random_index() % upper;
 }
 
-ScalarType uniform_random()
+inline ScalarType uniform_random()
 {
 #ifdef CUSTOM_UNIFORM_RANDOM_FUNCTION
 	return CUSTOM_UNIFORM_RANDOM_FUNCTION;
 #else
-	return rand()/static_cast<ScalarType>(RAND_MAX);
+	return std::rand()/((double)RAND_MAX+1);
 #endif
 }
 
-ScalarType gaussian_random()
+inline ScalarType gaussian_random()
 {
 #ifdef CUSTOM_GAUSSIAN_RANDOM_FUNCTION
 	return CUSTOM_GAUSSIAN_RANDOM_FUNCTION;
 #else
 	ScalarType x, y, radius;
 	do {
-		x = 2 * (rand() / ((double) RAND_MAX + 1)) - 1;
-		y = 2 * (rand() / ((double) RAND_MAX + 1)) - 1;
+		x = 2*(std::rand()/((double)RAND_MAX+1)) - 1;
+		y = 2*(std::rand()/((double)RAND_MAX+1)) - 1;
 		radius = (x * x) + (y * y);
 	} while ((radius >= 1.0) || (radius == 0.0));
-	radius = sqrt(-2 * log(radius) / radius);
+	radius = std::sqrt(-2 * std::log(radius) / radius);
 	x *= radius;
 	y *= radius;
 	return x;
@@ -54,9 +55,9 @@ ScalarType gaussian_random()
 }
 
 template <class RAI>
-void random_shuffle(RAI first, RAI last)
+inline void random_shuffle(RAI first, RAI last)
 {
-	std::random_shuffle(first,last);
+	std::random_shuffle(first,last,uniform_random_index_bounded);
 }
 
 }
diff --git a/src/shogun/lib/tapkee/methods.hpp b/src/shogun/lib/tapkee/methods.hpp
index 68ca1aae82a..00a93229ef5 100644
--- a/src/shogun/lib/tapkee/methods.hpp
+++ b/src/shogun/lib/tapkee/methods.hpp
@@ -467,7 +467,7 @@ class ImplementationBase
 	TapkeeOutput embedPassThru()
 	{
 		DenseMatrix feature_matrix =
-			dense_matrix_from_features(features, begin, end);
+			dense_matrix_from_features(features, current_dimension, begin, end);
 		return TapkeeOutput(feature_matrix.transpose(),tapkee::ProjectingFunction());
 	}
 
@@ -485,7 +485,7 @@ class ImplementationBase
 			               static_cast<ScalarType>((n_vectors-1)/3.0));
 
 		DenseMatrix data = 
-			dense_matrix_from_features(features, begin, end);
+			dense_matrix_from_features(features, current_dimension, begin, end);
 
 		DenseMatrix embedding(static_cast<IndexType>(target_dimension),n_vectors);
 		tsne::TSNE tsne;
@@ -501,7 +501,7 @@ class ImplementationBase
 			         static_cast<ScalarType>(1.0));
 
 		DenseMatrix embedding =
-			dense_matrix_from_features(features, begin, end);
+			dense_matrix_from_features(features, current_dimension, begin, end);
 
 		Neighbors neighbors = findNeighborsWith(plain_distance);
 
diff --git a/src/shogun/lib/tapkee/routines/eigendecomposition.hpp b/src/shogun/lib/tapkee/routines/eigendecomposition.hpp
index e2517bb39e2..536de5ee2e6 100644
--- a/src/shogun/lib/tapkee/routines/eigendecomposition.hpp
+++ b/src/shogun/lib/tapkee/routines/eigendecomposition.hpp
@@ -90,10 +90,9 @@ EigendecompositionResult eigendecomposition_impl_randomized(const MatrixType& wm
 	DenseMatrix O(wm.rows(), target_dimension+skip);
 	for (IndexType i=0; i<O.rows(); ++i)
 	{
-		for (IndexType j=0; j+1 < O.cols(); j++)
+		for (IndexType j=0; j<O.cols(); j++)
 		{
 			O(i,j) = tapkee::gaussian_random();
-			O(i,j+1) = tapkee::gaussian_random();
 		}
 	}
 	MatrixOperationType operation(wm);
diff --git a/src/shogun/lib/tapkee/routines/manifold_sculpting.hpp b/src/shogun/lib/tapkee/routines/manifold_sculpting.hpp
index 24ad5c18285..3d991350ffa 100644
--- a/src/shogun/lib/tapkee/routines/manifold_sculpting.hpp
+++ b/src/shogun/lib/tapkee/routines/manifold_sculpting.hpp
@@ -70,7 +70,8 @@ SparseMatrix neighbors_distances_matrix(RandomAccessIterator begin, RandomAccess
 {
 	const IndexType k = neighbors[0].size();
 	const IndexType n = neighbors.size();
-	assert((end-begin)==n);
+	if ((end-begin)!=n)
+		throw std::runtime_error("Wrong size");
 	SparseTriplets sparse_triplets;
 	sparse_triplets.reserve(k*n);
 	average_distance = 0;
@@ -92,7 +93,7 @@ SparseMatrix neighbors_distances_matrix(RandomAccessIterator begin, RandomAccess
 }
 
 SparseMatrixNeighborsPair angles_matrix_and_neighbors(const Neighbors& neighbors, 
-													const DenseMatrix& data)
+                                                      const DenseMatrix& data)
 {
 	const IndexType k = neighbors[0].size();
 	const IndexType n_vectors = data.cols();
@@ -158,7 +159,7 @@ ScalarType average_neighbor_distance(const DenseMatrix& data, const Neighbors& n
 }
 
 ScalarType compute_error_for_point(const IndexType index, const DenseMatrix& data,
-								   const DataForErrorFunc& error_func_data)
+                                   const DataForErrorFunc& error_func_data)
 {
 	IndexType k = error_func_data.distance_neighbors[0].size();
 	ScalarType error_value = 0;
@@ -213,10 +214,10 @@ ScalarType compute_error_for_point(const IndexType index, const DenseMatrix& dat
  * point
  */
 IndexType adjust_point_at_index(const IndexType index, DenseMatrix& data, 
-								const IndexType target_dimension, 
-								const ScalarType learning_rate,
-								const DataForErrorFunc& error_func_data,
-								ScalarType& point_error)
+                                const IndexType target_dimension, 
+                                const ScalarType learning_rate,
+                                const DataForErrorFunc& error_func_data,
+                                ScalarType& point_error)
 {
 	IndexType n_steps = 0;
 	ScalarType old_error, new_error;
diff --git a/src/shogun/lib/tapkee/utils/features.hpp b/src/shogun/lib/tapkee/utils/features.hpp
index 2d8d4a27908..3169f982d84 100644
--- a/src/shogun/lib/tapkee/utils/features.hpp
+++ b/src/shogun/lib/tapkee/utils/features.hpp
@@ -16,12 +16,13 @@ namespace tapkee_internal
 {
 
 template<class RandomAccessIterator, class FeaturesCallback>
-DenseMatrix dense_matrix_from_features(const FeaturesCallback& features,
-									const RandomAccessIterator& begin, 
-									const RandomAccessIterator& end)
+DenseMatrix dense_matrix_from_features(FeaturesCallback features,
+                                       IndexType dimension,
+                                       RandomAccessIterator begin, 
+                                       RandomAccessIterator end)
 {
-	DenseMatrix matrix(features.dimension(), end-begin);
-	DenseVector feature_vector(features.dimension());
+	DenseMatrix matrix(dimension, end-begin);
+	DenseVector feature_vector(dimension);
 
 	for (RandomAccessIterator iter=begin; iter!=end; ++iter)
 	{