From abb2ea110a237e2164ca92b32d2bf9c44607e0f8 Mon Sep 17 00:00:00 2001
From: Dmitry Vinokurov <dmitry.vinokurov@gmail.com>
Date: Fri, 20 Mar 2020 19:47:31 +0600
Subject: [PATCH 1/4] Issue #83 add suport of generic containers to VOI (#124)

* add utility classes for metaprogramming

* remove duplicate file

* enabling distance examples

* add support of generic containers to Tree class

* add support of generic containers to entropy and VOI

* add tests for entropy

* Update modules/distance/k-random/VOI.cpp

Co-Authored-By: jura-gresko <57671564+jura-gresko@users.noreply.github.com>

* some code cleaning

Co-authored-by: jura-gresko <57671564+jura-gresko@users.noreply.github.com>
---
 examples/CMakeLists.txt                       |   2 +-
 examples/distance_examples/CMakeLists.txt     |   8 +-
 .../distance_examples/entropy_example.cpp     |   2 +-
 modules/distance/k-random/VOI.cpp             | 154 ++++++++--------
 modules/distance/k-random/VOI.hpp             |  77 ++++----
 modules/space/tree.cpp                        |   3 +-
 modules/space/tree.hpp                        |   3 +-
 modules/utils/type_traits.hpp                 | 157 +++++++++++++++++
 tests/CMakeLists.txt                          |   3 +-
 tests/distance_tests/CMakeLists.txt           |  31 ++++
 tests/distance_tests/voi_tests.cpp            | 164 ++++++++++++++++++
 11 files changed, 490 insertions(+), 114 deletions(-)
 create mode 100644 modules/utils/type_traits.hpp
 create mode 100644 tests/distance_tests/CMakeLists.txt
 create mode 100644 tests/distance_tests/voi_tests.cpp

diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt
index 899f89b8..57dfc518 100644
--- a/examples/CMakeLists.txt
+++ b/examples/CMakeLists.txt
@@ -1,6 +1,6 @@
 add_subdirectory(correlation_examples)
 add_subdirectory(energies_examples)
-#add_subdirectory(distance_examples)
+add_subdirectory(distance_examples)
 add_subdirectory(dnn_examples)
 add_subdirectory(ensemble_examples)
 #add_subdirectory(graph_pyhton_assets)
diff --git a/examples/distance_examples/CMakeLists.txt b/examples/distance_examples/CMakeLists.txt
index 142feee5..59115998 100644
--- a/examples/distance_examples/CMakeLists.txt
+++ b/examples/distance_examples/CMakeLists.txt
@@ -19,7 +19,7 @@ include_directories(${Boost_INCLUDE_DIRS})
 file(COPY "assets" DESTINATION ".")
 
 file(GLOB EXAMPLE_SRCS RELATIVE ${CMAKE_CURRENT_SOURCE_DIR} *.cpp) 
-
+find_package(JPEG) 
 #Run through each source
 foreach(exampleSrc ${EXAMPLE_SRCS})
 
@@ -47,8 +47,10 @@ foreach(exampleSrc ${EXAMPLE_SRCS})
 endforeach(exampleSrc)
 
 if(UNIX)
-	target_link_libraries(earth_mover_distance_example /usr/local/lib/libjpeg.a) 
-	target_link_libraries(earth_mover_distance_2_example /usr/local/lib/libjpeg.a) 
+	target_link_libraries(earth_mover_distance_example JPEG::JPEG)
+	target_link_libraries(earth_mover_distance_2_example JPEG::JPEG)
+#	target_link_libraries(earth_mover_distance_example /usr/local/lib/libjpeg.a) 
+#	target_link_libraries(earth_mover_distance_2_example /usr/local/lib/libjpeg.a) 
 endif(UNIX)
 
 if(CMAKE_SYSTEM_NAME MATCHES Windows)
diff --git a/examples/distance_examples/entropy_example.cpp b/examples/distance_examples/entropy_example.cpp
index ff79ab8f..dac3d260 100644
--- a/examples/distance_examples/entropy_example.cpp
+++ b/examples/distance_examples/entropy_example.cpp
@@ -88,7 +88,7 @@ int main() {
 	std::cout << "Variation of Information, normalized Variation of Information:" << std::endl;
 	
 	std::cout << "VOI = " << metric::variationOfInformation(v1, v2) << std::endl;
-	std::cout << "VOI (Manhatten) = " << metric::variationOfInformation<double, metric::Manhatten<double>>(v1, v2) << std::endl;
+	std::cout << "VOI (Manhatten) = " << metric::variationOfInformation<std::vector<std::vector<double>>, metric::Manhatten<double>>(v1, v2) << std::endl;
 	std::cout << "VOI norm = " << metric::variationOfInformation_normalized(v1, v2) << std::endl;
 
 	// functor
diff --git a/modules/distance/k-random/VOI.cpp b/modules/distance/k-random/VOI.cpp
index eb56eafc..d096a886 100644
--- a/modules/distance/k-random/VOI.cpp
+++ b/modules/distance/k-random/VOI.cpp
@@ -13,6 +13,7 @@ Copyright (c) 2019 Panda Team
 #include <iostream>
 #include <random>
 #include <unordered_set>
+#include <utility>
 #include <vector>
 
 #include <boost/functional/hash.hpp>
@@ -27,9 +28,10 @@ Copyright (c) 2019 Panda Team
 namespace metric {
 
 namespace {
-    template <typename T>
-    void add_noise(std::vector<std::vector<T>>& data)
+    template <typename C>
+    void add_noise(C & data)
     {
+        using T = type_traits::underlaying_type_t<C>;
         std::random_device rd;
         std::mt19937 gen(rd());
         std::normal_distribution<T> dis(0, 1);
@@ -69,14 +71,13 @@ namespace {
         }
         std::cout << "]" << std::endl;
     }
-    template <typename T>
-    void combine(
-        const std::vector<std::vector<T>>& X, const std::vector<std::vector<T>>& Y, std::vector<std::vector<T>>& XY)
+    template <typename C1, typename C2, typename T=type_traits::underlaying_type_t<C1>>
+    std::vector<std::vector<T>> combine(const C1& X, const C2& Y)
     {
         std::size_t N = X.size();
         std::size_t dx = X[0].size();
         std::size_t dy = Y[0].size();
-        XY.resize(N);
+        std::vector<std::vector<T>> XY(N);
         for (std::size_t i = 0; i < N; i++) {
             XY[i].resize(dx + dy);
             std::size_t k = 0;
@@ -87,9 +88,10 @@ namespace {
                 XY[i][k] = Y[i][j];
             }
         }
+        return XY;
     }
-    template <typename T>
-    std::vector<T> unique(const std::vector<T>& data)
+    template <typename Container, typename T = type_traits::underlaying_type_t<Container>>
+    std::vector<T> unique(const Container& data)
     {
         std::unordered_set<std::size_t> hashes;
         std::vector<T> result;
@@ -104,11 +106,11 @@ namespace {
 
 template <typename Container, typename Metric, typename L>
 double entropy(
-    std::vector<Container> data, std::size_t k, L logbase, Metric metric)
+    const Container & data, std::size_t k, L logbase, const Metric & metric)
 {
-    using T = typename Container::value_type;
-
-    if (data.empty() || data[0].empty()) {
+    using T_underlaying = type_traits::underlaying_type_t<Container>;
+    using T = type_traits::index_value_type_t<Container>; 
+    if (data.size() == 0 || data[0].size() == 0) {
         return 0;
     }
     if (data.size() < k + 1)
@@ -120,17 +122,17 @@ double entropy(
     double two = 2.0;  // this is in order to make types match the log template function
     double cb = d * log(logbase, two);
 
-    if constexpr (!std::is_same<Metric, typename metric::Chebyshev<T>>::value) {
-        if constexpr (std::is_same<Metric, typename metric::Euclidian<T>>::value) {
+    if constexpr (!std::is_same<Metric, typename metric::Chebyshev<T_underlaying>>::value) {
+        if constexpr (std::is_same<Metric, typename metric::Euclidian<T_underlaying>>::value) {
             p = 2;
-        } else if constexpr (std::is_same<Metric, typename metric::P_norm<T>>::value) {
+        } else if constexpr (std::is_same<Metric, typename metric::P_norm<T_underlaying>>::value) {
             p = metric.p;
         }
         cb = cb + d * log(logbase, std::tgamma(1 + 1 / p)) - log(logbase, std::tgamma(1 + d / p));
     }
 
     //add_noise(data); // TODO test
-    metric::Tree<Container, Metric> tree(data, -1, metric);
+    metric::Tree<T, Metric> tree(data, -1, metric);
     double entropyEstimate = boost::math::digamma(N) - boost::math::digamma(k) + cb + d * log(logbase, two);
     for (std::size_t i = 0; i < N; i++) {
         auto res = tree.knn(data[i], k + 1);
@@ -142,10 +144,13 @@ double entropy(
 // Kozachenko-Leonenko estimator based on https://hal.archives-ouvertes.fr/hal-00331300/document (Shannon diff. entropy,
 // q = 1)
 
-template <typename T, typename Metric = metric::Euclidian<T>, typename L = T>  // TODO check if L = T is correct
-typename std::enable_if<!std::is_integral<T>::value, T>::type entropy_kl(
-    std::vector<std::vector<T>> data, std::size_t k = 3, L logbase = 2, Metric metric = Metric())
+template <typename Container, typename Metric = metric::Euclidian<type_traits::underlaying_type_t<Container>>,
+    typename L = type_traits::underlaying_type_t<Container>>  // TODO check if L = T is correct
+std::enable_if_t<!type_traits::is_container_of_integrals_v<Container>, L>
+entropy_kl( const Container & data, std::size_t k = 3, L logbase = 2, const Metric &metric = Metric())
 {
+    using T1 = type_traits::index_value_type_t<Container>;
+    using T = type_traits::underlaying_type_t<Container>;
     if (data.empty() || data[0].empty())
         return 0;
     if (data.size() < k + 1)
@@ -153,7 +158,7 @@ typename std::enable_if<!std::is_integral<T>::value, T>::type entropy_kl(
     if constexpr (!std::is_same<Metric, typename metric::Euclidian<T>>::value)
         throw std::logic_error("entropy function is now implemented only for Euclidean distance");
 
-    metric::Tree<std::vector<T>, Metric> tree(data, -1, metric);
+    metric::Tree<T1, Metric> tree(data, -1, metric);
 
     size_t N = data.size();
     size_t m = data[0].size();
@@ -172,8 +177,8 @@ typename std::enable_if<!std::is_integral<T>::value, T>::type entropy_kl(
     return sum;
 }
 
-template <typename T>
-std::pair<std::vector<double>, std::vector<std::vector<T>>> pluginEstimator(const std::vector<std::vector<T>>& Y)
+template <typename Container, typename T>
+std::pair<std::vector<double>, Container> pluginEstimator(const Container & Y)
 {
     std::vector<std::vector<T>> uniqueVal = unique(Y);
     std::vector<double> counts(uniqueVal.size());
@@ -189,10 +194,14 @@ std::pair<std::vector<double>, std::vector<std::vector<T>>> pluginEstimator(cons
     return std::make_pair(counts, uniqueVal);
 }
 
-template <typename T, typename Metric>
-typename std::enable_if<!std::is_integral<T>::value, T>::type mutualInformation(
-    const std::vector<std::vector<T>>& Xc, const std::vector<std::vector<T>>& Yc, int k, Metric metric, int version)
+template <typename Container, typename Metric>
+typename std::enable_if_t<!type_traits::is_container_of_integrals_v<Container>,
+                          type_traits::underlaying_type_t<Container>>
+mutualInformation(
+                  const Container & Xc, const Container & Yc, int k, const Metric & metric, int version)
 {
+    using T = type_traits::underlaying_type_t<Container>;
+    using V = type_traits::index_value_type_t<Container>;
     T N = Xc.size();
 
     if (N < k + 1 || Yc.size() < k + 1)
@@ -202,15 +211,14 @@ typename std::enable_if<!std::is_integral<T>::value, T>::type mutualInformation(
     auto Y = Yc;
     add_noise(X);
     add_noise(Y);
-    std::vector<std::vector<T>> XY;
-    combine(X, Y, XY);
+    std::vector<std::vector<T>> XY = combine(X, Y);
     metric::Tree<std::vector<T>, Metric> tree(XY, -1, metric);
     auto entropyEstimate = boost::math::digamma(k) + boost::math::digamma(N);
     if (version == 2) {
         entropyEstimate -= 1 / static_cast<double>(k);
     }
 
-    metric::Tree<std::vector<T>, Metric> xTree(X, -1, metric);
+    metric::Tree<V, Metric> xTree(X, -1, metric);
 
     for (std::size_t i = 0; i < N; i++) {
         auto res = tree.knn(XY[i], k + 1);
@@ -239,83 +247,87 @@ typename std::enable_if<!std::is_integral<T>::value, T>::type mutualInformation(
     return entropyEstimate;
 }
 
-template <typename T>
-typename std::enable_if<std::is_integral<T>::value, T>::type mutualInformation(
-    const std::vector<std::vector<T>>& Xc, const std::vector<std::vector<T>>& Yc, T logbase)
+template <typename Container, typename T>
+std::enable_if_t<type_traits::is_container_of_integrals_v<Container>, T>
+mutualInformation(const Container& Xc, const Container& Yc, T logbase)
 {
-    std::vector<std::vector<T>> XY;
-    combine(Xc, Yc, XY);
-    return entropy<T>(Xc, logbase)
-        + entropy<T>(Yc,
+    std::vector<std::vector<T>> XY = combine(Xc, Yc);
+    return entropy(Xc, logbase)
+        + entropy(Yc,
             logbase)  // entropy overload for integers is not implemented yet
-        - entropy<T>(XY, logbase);
+        - entropy(XY, logbase);
 }
 
-template <typename T, typename Metric>
-typename std::enable_if<!std::is_integral<T>::value, T>::type variationOfInformation(
-    const std::vector<std::vector<T>>& Xc, const std::vector<std::vector<T>>& Yc, int k, T logbase)
+template <typename C, typename Metric, typename T>
+typename std::enable_if_t<!type_traits::is_container_of_integrals_v<C>, T>
+variationOfInformation(const C& Xc, const C& Yc, int k, T logbase)
 {
-    return entropy<std::vector<T>, Metric>(Xc, k, logbase, Metric()) + entropy<std::vector<T>, Metric>(Yc, k, logbase, Metric())
-        - 2 * mutualInformation<T>(Xc, Yc, k);
+    return entropy<C, Metric>(Xc, k, logbase, Metric()) + entropy<C, Metric>(Yc, k, logbase, Metric())
+        - 2 * mutualInformation<C>(Xc, Yc, k);
 }
 
-template <typename T>
-typename std::enable_if<!std::is_integral<T>::value, T>::type variationOfInformation_normalized(
-    const std::vector<std::vector<T>>& Xc, const std::vector<std::vector<T>>& Yc, int k, T logbase)
+template <typename C, typename T>
+typename std::enable_if_t<!type_traits::is_container_of_integrals_v<C>, T>
+variationOfInformation_normalized(
+        const C& Xc, const C& Yc, int k, T logbase)
 {
     using Cheb = metric::Chebyshev<T>;
-    auto mi = mutualInformation<T>(Xc, Yc, k);
-    return 1 - (mi / (entropy<std::vector<T>, Cheb>(Xc, k, logbase, Cheb()) + entropy<std::vector<T>, Cheb>(Yc, k, logbase, Cheb()) - mi));
+    auto mi = mutualInformation<C>(Xc, Yc, k);
+    return 1 - (mi / (entropy(Xc, k, logbase, Cheb()) + entropy(Yc, k, logbase, Cheb()) - mi));
+}
+
+template <typename V>
+template <typename C, typename El>
+typename std::enable_if_t<!std::is_integral_v<El>, V>  // only real values are accepted
+VOI<V>::operator()(const C& a, const C& b) const
+{
+    using Cheb = metric::Chebyshev<El>;
+    return entropy(a, k, logbase, Cheb()) + entropy(b, k, logbase, Cheb())
+        - 2 * mutualInformation<C>(a, b, k);
 }
 
 template <typename V>
-template <template <class, class> class Container, class Allocator_inner, class Allocator_outer, class El>
-typename std::enable_if<!std::is_integral<El>::value, V>::type VOI<V>::operator()(
-    const Container<Container<El, Allocator_inner>, Allocator_outer>& a,
-    const Container<Container<El, Allocator_inner>, Allocator_outer>& b) const
+template <typename C, typename El>
+typename std::enable_if_t<std::is_integral_v<El>, V>  // only real values are accepted
+VOI<V>::operator()(const C& a, const C& b) const 
 {
     using Cheb = metric::Chebyshev<El>;
-    return entropy<std::vector<El>, Cheb>(a, k, logbase, Cheb()) + entropy<std::vector<El>, Cheb>(b, k, logbase, Cheb())
-        - 2 * mutualInformation<El>(a, b, k);
+    return entropy(a, k, logbase, Cheb()) + entropy(b, k, logbase, Cheb())
+        - 2 * mutualInformation<C>(a, b, k);
 }
 
 template <typename V>
-template <template <class, class> class Container, class Allocator_inner, class Allocator_outer, class El>
-typename std::enable_if<!std::is_integral<El>::value, V>::type VOI_normalized<V>::operator()(
-    const Container<Container<El, Allocator_inner>, Allocator_outer>& a,
-    const Container<Container<El, Allocator_inner>, Allocator_outer>& b) const
+template <typename C, typename El>
+typename std::enable_if_t<!std::is_integral_v<El>, V>  // only real values are accepted
+VOI_normalized<V>::operator()(const C& a, const C& b) const
 {
     using Cheb = metric::Chebyshev<El>;
-    auto mi = mutualInformation<El>(a, b, this->k);
+    auto mi = mutualInformation(a, b, this->k);
     return 1
         - (mi
-            / (entropy<std::vector<El>, Cheb>(a, this->k, this->logbase, Cheb())
-                + entropy<std::vector<El>, Cheb>(b, this->k, this->logbase, Cheb()) - mi));
+            / (entropy(a, this->k, this->logbase, Cheb())
+                + entropy(b, this->k, this->logbase, Cheb()) - mi));
 }
 
 // VOI based on Kozachenko-Leonenko entropy estimator
 
 template <typename V>
-template <template <class, class> class Container, class Allocator_inner, class Allocator_outer, class El>
-typename std::enable_if<!std::is_integral<El>::value, V>::type VOI_kl<V>::operator()(
-    const Container<Container<El, Allocator_inner>, Allocator_outer>& a,
-    const Container<Container<El, Allocator_inner>, Allocator_outer>& b) const
+template <typename C, typename El>
+typename std::enable_if_t<!std::is_integral_v<El>, V>  // only real values are accepted
+    VOI_kl<V>::operator()(const C& a, const C& b) const
 {
-    std::vector<std::vector<El>> ab;
-    combine(a, b, ab);
+    std::vector<std::vector<El>> ab = combine(a, b);
     return 2 * entropy_kl(ab) - entropy_kl(a) - entropy_kl(b);
 }
 
 template <typename V>
-template <template <class, class> class Container, class Allocator_inner, class Allocator_outer, class El>
-typename std::enable_if<!std::is_integral<El>::value, V>::type VOI_normalized_kl<V>::operator()(
-    const Container<Container<El, Allocator_inner>, Allocator_outer>& a,
-    const Container<Container<El, Allocator_inner>, Allocator_outer>& b) const
+template <typename C, typename El>
+typename std::enable_if_t<!std::is_integral_v<El>, V>  // only real values are accepted
+VOI_normalized_kl<V>::operator()(const C& a, const C& b) const
 {
     auto entropy_a = entropy_kl(a);
     auto entropy_b = entropy_kl(b);
-    std::vector<std::vector<El>> ab;
-    combine(a, b, ab);
+    std::vector<std::vector<El>> ab = combine(a, b);
     auto joint_entropy = entropy_kl(ab);
     auto mi = entropy_a + entropy_b - joint_entropy;
     return 1 - (mi / (entropy_a + entropy_b - mi));
diff --git a/modules/distance/k-random/VOI.hpp b/modules/distance/k-random/VOI.hpp
index a31a3458..0fecaf52 100644
--- a/modules/distance/k-random/VOI.hpp
+++ b/modules/distance/k-random/VOI.hpp
@@ -12,9 +12,8 @@ Copyright (c) 2019 Panda Team
 #include "../k-related/Standards.hpp"
 #include <type_traits>
 #include <vector>
-
+#include "../../../modules/utils/type_traits.hpp"
 namespace metric {
-
 /**
  * @brief Continuous entropy estimator
  *
@@ -24,9 +23,11 @@ namespace metric {
  * @param metric
  * @return value of entropy estimation of the data 
  */
-template <typename Container, typename Metric = metric::Euclidian<typename Container::value_type>, typename L = double>
-double entropy(
-    std::vector<Container> data, std::size_t k = 3, L logbase = 2, Metric metric = Metric());
+template<typename Container,
+         typename Metric = metric::Euclidian<type_traits::underlaying_type_t<Container>>,
+         typename L = double
+         >
+double entropy(const Container & data, std::size_t k = 3, L logbase = 2, const Metric & metric = Metric());
 
 /**
  * @brief
@@ -34,8 +35,8 @@ double entropy(
  * @param Y
  * @return
  */
-template <typename T>
-std::pair<std::vector<double>, std::vector<std::vector<T>>> pluginEstimator(const std::vector<std::vector<T>>& Y);
+template <typename Container, typename T = type_traits::underlaying_type_t<Container>>
+std::pair<std::vector<double>, std::vector<std::vector<T>>> pluginEstimator(const Container& Y);
 
 /**
  * @brief
@@ -47,13 +48,15 @@ std::pair<std::vector<double>, std::vector<std::vector<T>>> pluginEstimator(cons
  * @param version
  * @return
  */
-template <typename T, typename Metric = metric::Chebyshev<T>>
-typename std::enable_if<!std::is_integral<T>::value, T>::type mutualInformation(const std::vector<std::vector<T>>& Xc,
-    const std::vector<std::vector<T>>& Yc, int k = 3, Metric metric = Metric(), int version = 2);
+template <typename Container, typename Metric = metric::Chebyshev<type_traits::underlaying_type_t<Container>>>
+typename std::enable_if_t<!type_traits::is_container_of_integrals_v<Container>,
+                          type_traits::underlaying_type_t<Container>>
+mutualInformation(const Container& Xc,
+    const Container & Yc, int k = 3, const Metric & metric = Metric(), int version = 2);
 
-template <typename T>
-typename std::enable_if<std::is_integral<T>::value, T>::type mutualInformation(
-    const std::vector<std::vector<T>>& Xc, const std::vector<std::vector<T>>& Yc, T logbase = 2.0);
+template <typename Container, typename T = type_traits::underlaying_type_t<Container> >
+std::enable_if_t<type_traits::is_container_of_integrals_v<Container>, T>
+mutualInformation(const Container& Xc, const Container& Yc, T logbase = 2.0);
 
 /**
  * @brief
@@ -64,9 +67,11 @@ typename std::enable_if<std::is_integral<T>::value, T>::type mutualInformation(
  * @param logbase
  * @return
  */
-template <typename T = double, typename Metric = metric::Chebyshev<T>>
-typename std::enable_if<!std::is_integral<T>::value, T>::type variationOfInformation(
-    const std::vector<std::vector<T>>& Xc, const std::vector<std::vector<T>>& Yc, int k = 3, T logbase = 2.0);
+template <typename C, typename Metric = metric::Chebyshev<type_traits::underlaying_type_t<C>>,
+          typename T = type_traits::underlaying_type_t<C>>
+typename std::enable_if_t<!type_traits::is_container_of_integrals_v<C>, T>
+variationOfInformation(const C& Xc, const C& Yc, int k = 3, T logbase = 2.0);
+
 
 /**
  * @brief
@@ -77,9 +82,10 @@ typename std::enable_if<!std::is_integral<T>::value, T>::type variationOfInforma
  * @param logbase
  * @return
  */
-template <typename T = double>
-typename std::enable_if<!std::is_integral<T>::value, T>::type variationOfInformation_normalized(
-    const std::vector<std::vector<T>>& Xc, const std::vector<std::vector<T>>& Yc, int k = 3, T logbase = 2.0);
+    template <typename C, typename T = type_traits::underlaying_type_t<C>>
+    typename std::enable_if_t<!type_traits::is_container_of_integrals_v<C>,T>
+    variationOfInformation_normalized(
+    const C & Xc, const C& Yc, int k = 3, T logbase = 2.0);
 
 /**
  * @brief
@@ -112,10 +118,14 @@ struct VOI {
      * @param b second container
      * @return variation of information between a and b
      */
-    template <template <class, class> class Container, class Allocator_inner, class Allocator_outer, class El>
-    typename std::enable_if<!std::is_integral<El>::value, V>::type  // only real values are accepted
-    operator()(const Container<Container<El, Allocator_inner>, Allocator_outer>& a,
-        const Container<Container<El, Allocator_inner>, Allocator_outer>& b) const;
+
+    template <typename Container, typename T = type_traits::underlaying_type_t<Container>>
+    typename std::enable_if_t<!std::is_integral_v<T>, V>  // only real values are accepted
+    operator()(const Container& a, const Container& b) const;
+
+    template <typename Container, typename T = type_traits::underlaying_type_t<Container>>
+    typename std::enable_if_t<std::is_integral_v<T>, V>  // only real values are accepted
+    operator()(const Container& a, const Container& b) const;
 
     // TODO add support of 1D random values passed in simple containers
 };
@@ -149,10 +159,9 @@ struct VOI_normalized : VOI<V> {
      * @param b second container
      * @return varition of information between a and b
      */
-    template <template <class, class> class Container, class Allocator_inner, class Allocator_outer, class El>
-    typename std::enable_if<!std::is_integral<El>::value, V>::type  // only real values are accepted
-    operator()(const Container<Container<El, Allocator_inner>, Allocator_outer>& a,
-        const Container<Container<El, Allocator_inner>, Allocator_outer>& b) const;
+    template <typename Container, typename T = type_traits::underlaying_type_t<Container>>
+    typename std::enable_if_t<!std::is_integral_v<T>, V>  // only real values are accepted
+    operator()(const Container& a, const Container& b) const;
 
     // TODO add support of 1D random values passed in simple containers
 };
@@ -192,10 +201,9 @@ struct VOI_kl {
      * @param b second container
      * @return variation of information
      */
-    template <template <class, class> class Container, class Allocator_inner, class Allocator_outer, class El>
-    typename std::enable_if<!std::is_integral<El>::value, V>::type  // only real values are accepted
-    operator()(const Container<Container<El, Allocator_inner>, Allocator_outer>& a,
-        const Container<Container<El, Allocator_inner>, Allocator_outer>& b) const;
+    template <typename Container, typename T = type_traits::underlaying_type_t<Container>>
+    typename std::enable_if_t<!std::is_integral_v<T>, V>  // only real values are accepted
+    operator()(const Container& a, const Container& b) const;
 
     // TODO add support of 1D random values passed in simple containers
 };
@@ -228,10 +236,9 @@ struct VOI_normalized_kl : VOI_kl<V> {
      * @param b second value
      * @return normalized value of the variation of information
      */
-    template <template <class, class> class Container, class Allocator_inner, class Allocator_outer, class El>
-    typename std::enable_if<!std::is_integral<El>::value, V>::type  // only real values are accepted
-    operator()(const Container<Container<El, Allocator_inner>, Allocator_outer>& a,
-        const Container<Container<El, Allocator_inner>, Allocator_outer>& b) const;
+    template <typename Container, typename T = type_traits::underlaying_type_t<Container>>
+    typename std::enable_if_t<!std::is_integral_v<T>, V>  // only real values are accepted
+    operator()(const Container& a, const Container& b) const;
 
     // TODO add support of 1D random values passed in simple containers
 };
diff --git a/modules/space/tree.cpp b/modules/space/tree.cpp
index f2ac5e48..04b65af6 100644
--- a/modules/space/tree.cpp
+++ b/modules/space/tree.cpp
@@ -324,7 +324,8 @@ Tree<recType, Metric>::Tree(const recType& p, int truncateArg /*=-1*/, Metric d)
 
 /*** constructor: with a vector data records **/
 template <class recType, class Metric>
-Tree<recType, Metric>::Tree(const std::vector<recType>& p, int truncateArg /*=-1*/, Metric d)
+template<typename C>
+Tree<recType, Metric>::Tree(const C& p, int truncateArg /*=-1*/, Metric d)
     : metric_(d)
 {
     min_scale = 1000;
diff --git a/modules/space/tree.hpp b/modules/space/tree.hpp
index 9d356883..7026588e 100644
--- a/modules/space/tree.hpp
+++ b/modules/space/tree.hpp
@@ -149,7 +149,8 @@ class Tree {
      * @param truncate truncate paramter
      * @param d metric object
      */
-    Tree(const std::vector<recType>& p, int truncate = -1, Metric d = Metric());  // with a vector of data records
+    template<typename Container>
+    Tree(const Container& p, int truncate = -1, Metric d = Metric());  // with a vector of data records
 
     /**
      * @brief Destroy the Tree object
diff --git a/modules/utils/type_traits.hpp b/modules/utils/type_traits.hpp
new file mode 100644
index 00000000..04819199
--- /dev/null
+++ b/modules/utils/type_traits.hpp
@@ -0,0 +1,157 @@
+/*
+  This Source Code Form is subject to the terms of the Mozilla Public
+  License, v. 2.0. If a copy of the MPL was not distributed with this
+  file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+  Copyright (c) 2020 Panda Team
+*/
+
+#ifndef METRIC_UTILS_TYPE_TRAITS_HPP
+#define METRIC_UTILS_TYPE_TRAITS_HPP
+#include <utility>
+#include <type_traits>
+
+namespace metric::type_traits {
+/**
+ * @brief if T is a container and implemented operator[] value is true, otherwise value is false
+ *
+ * @tparam T checking type
+ *
+ */
+template <typename T>
+class is_has_index_operator {
+    struct nil_t {
+    };
+    template <typename U>
+    static constexpr auto test(U*) -> decltype(std::declval<U&>()[0]);
+    template <typename>
+    static constexpr auto test(...) -> nil_t;
+
+public:
+    using type = typename std::decay<decltype(test<T>(nullptr))>::type;
+    static const bool value = !std::is_same<type, nil_t>::value;
+};
+
+/**
+ * @brief if T is a class that has resize() method value is true, otherwise value is false
+ *
+ * @tparam T checking type
+ */
+template <typename T>
+class is_has_resize_method {
+    struct nil_t {
+    };
+    template <typename U>
+    static constexpr auto test(U*) -> decltype(std::declval<U&>().resize(0));
+    template <typename>
+    static constexpr auto test(...) -> nil_t;
+
+public:
+    using type = typename std::decay<decltype(test<T>(nullptr))>::type;
+    static const bool value = !std::is_same<type, nil_t>::value;
+};
+
+/**
+ * @brief extract returning type of operator[] in container, in case of STL containers is equivalent ot value_type
+ * for example: underlaying_type<std::vector<int>> == int,
+ *              underlaying_type<std::vector<std::vector<int>>> == std::vector<int>,
+ 
+ * @tparam T Container type
+ */
+template <typename T>
+using index_value_type_t = typename is_has_index_operator<T>::type;
+
+/**
+ * @brief helper for is_has_index_operator metafunction
+ */
+template <typename T>
+inline constexpr bool is_has_index_operator_v = is_has_index_operator<T>::value;
+
+/**
+ * @brief helper for is_has_resize_method metafunction
+ */
+template <typename T>
+inline constexpr bool is_has_resize_method_v = is_has_resize_method<T>::value;
+
+/**
+ * @brief if T is a container supported [] operator and underlaying type of container is integral type
+ * value is true, otherwise value is false
+ */
+template <typename T>
+struct is_container_of_integrals {
+    static constexpr auto test()
+    {
+        if constexpr (is_has_index_operator_v<T>) {
+            using T1 = index_value_type_t<T>;
+            if constexpr (is_has_index_operator_v<T1>) {
+                using T2 = index_value_type_t<T1>;
+                return std::is_integral_v<T2>;
+            } else {
+                return std::is_integral_v<T1>;
+            }
+        } else {
+            return false;
+        }
+    }
+    static const bool value = test();
+};
+
+/**
+ * @brief helper for is_container_of_integrals metafunction
+ */
+template <typename T>
+inline constexpr bool is_container_of_integrals_v = is_container_of_integrals<T>::value;
+
+template <typename T>
+struct underlaying_type_impl0 {
+    using type = typename std::decay<T>::type;
+};
+template <typename T>
+using underlaying0_t = typename underlaying_type_impl0<T>::type;
+
+template <typename T>
+struct underlaying_type_impl1 {
+    using type = typename std::decay<typename is_has_index_operator<T>::type>::type;
+};
+template <typename T>
+using underlaying1_t = typename underlaying_type_impl1<T>::type;
+
+template <typename T>
+struct underlaying_type_impl2 {
+    using type =
+        typename std::decay<typename is_has_index_operator<typename is_has_index_operator<T>::type>::type>::type;
+};
+template <typename T>
+using underlaying2_t = typename underlaying_type_impl2<T>::type;
+
+/**
+ * @brief extract underlaying type of container,
+ * for example: underlaying_type<std::vector<int>> == int,
+ *              underlaying_type<std::vector<std::vector<int>>> == int,
+ *
+ */
+template <typename T>
+struct underlaying_type {
+    static constexpr auto level() -> int
+    {
+        if constexpr (is_has_index_operator_v<T>) {
+            using T1 = index_value_type_t<T>;
+            if constexpr (is_has_index_operator_v<T1>) {
+                return 2;
+            }
+            return 1;
+        }
+        return 0;
+    }
+    using type = std::conditional_t<level() == 2, underlaying2_t<T>,
+        std::conditional_t<level() == 1, underlaying1_t<T>, underlaying0_t<T>>>;
+};
+
+/**
+ * @brief helper for underlaying_type metafunction
+ */
+template <typename T>
+using underlaying_type_t = typename underlaying_type<T>::type;
+
+};
+#endif
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index ac100bc6..761f81e4 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -6,4 +6,5 @@ add_subdirectory(ensembles_tests)
 add_subdirectory(mapping_tests)
 add_subdirectory(space_tests)
 add_subdirectory(sparsification_tests)
-add_subdirectory(utils)
\ No newline at end of file
+add_subdirectory(utils)
+add_subdirectory(distance_tests)
\ No newline at end of file
diff --git a/tests/distance_tests/CMakeLists.txt b/tests/distance_tests/CMakeLists.txt
new file mode 100644
index 00000000..a9904423
--- /dev/null
+++ b/tests/distance_tests/CMakeLists.txt
@@ -0,0 +1,31 @@
+set(BOOST_TEST_LIB ${Boost_LIBRARIES})
+
+if(UNIX)
+  set(THREADS_PREFER_PTHREAD_FLAG ON)
+  find_package(Threads REQUIRED)
+endif(UNIX)
+
+AUX_SOURCE_DIRECTORY(. TEST_SRCS)
+#Run through each source
+foreach(testSrc ${TEST_SRCS})
+
+  #Extract the filename without an extension (NAME_WE)
+  get_filename_component(testName ${testSrc} NAME_WE)
+
+  #Add compile target
+  add_executable(${testName} ${testSrc})
+
+  #link to Boost libraries AND your targets and dependencies
+  if(WIN32)
+    target_compile_options(${testName} PRIVATE /D BOOST_ALL_NO_LIB=1)
+  endif()
+  target_link_libraries(${testName} ${BOOST_TEST_LIB})
+
+  if(UNIX)
+    target_link_libraries(${testName} Threads::Threads)
+# uncomment next two lines if want to enable address sanitizer in debug mode
+#    set (CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -fno-omit-frame-pointer -fsanitize=address")
+#    set (CMAKE_LINKER_FLAGS_DEBUG "${CMAKE_LINKER_FLAGS_DEBUG} -fno-omit-frame-pointer -fsanitize=address")
+  endif(UNIX)
+  add_test(NAME ${testName}  COMMAND ${testName} )
+endforeach(testSrc)
diff --git a/tests/distance_tests/voi_tests.cpp b/tests/distance_tests/voi_tests.cpp
new file mode 100644
index 00000000..e78bd8c8
--- /dev/null
+++ b/tests/distance_tests/voi_tests.cpp
@@ -0,0 +1,164 @@
+/*
+This Source Code Form is subject to the terms of the Mozilla Public
+License, v. 2.0. If a copy of the MPL was not distributed with this
+file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+Copyright (c) 2020 Panda Team
+*/
+
+#include "modules/utils/type_traits.hpp"
+#include <stdexcept>
+#define BOOST_TEST_DYN_LINK
+#define BOOST_TEST_MODULE voi_test
+#include <boost/test/unit_test.hpp>
+
+#include "modules/distance/k-random/VOI.hpp"
+#include "modules/distance/k-structured/Edit.hpp"
+#include "modules/distance/k-related/Standards.hpp"
+#include <vector>
+#include <deque>
+#include <array>
+#include <numeric>
+
+template <typename T>
+bool float_eq(T l, T r) {
+    return std::abs(l -r) <= std::numeric_limits<T>::epsilon();
+}
+BOOST_AUTO_TEST_CASE(entropy) {
+    std::vector<std::vector<double>>  v1 = { { 5, 5 }, { 2, 2 }, { 3, 3 }, { 5, 1 } };
+    std::vector<std::deque<double>>   v2 = { { 5, 5 }, { 2, 2 }, { 3, 3 }, { 5, 1 } };
+    std::vector<std::array<double,2>> v3 = { { 5, 5 }, { 2, 2 }, { 3, 3 }, { 5, 1 } };
+    std::deque<std::vector<double>>   v4 = { { 5, 5 }, { 2, 2 }, { 3, 3 }, { 5, 1 } };
+    std::deque<std::array<double, 2>> v7 = { { 5, 5 }, { 2, 2 }, { 3, 3 }, { 5, 1 } };
+    BOOST_TEST(float_eq(metric::entropy(v1, 3, 2.0, metric::Chebyshev<double>()), 7.625814583693912));
+    BOOST_TEST(float_eq(metric::entropy(v2, 3, 2.0, metric::Chebyshev<double>()), 7.625814583693912));
+    BOOST_TEST(float_eq(metric::entropy(v3, 3, 2.0, metric::Chebyshev<double>()), 7.625814583693912));
+    BOOST_TEST(float_eq(metric::entropy(v4, 3, 2.0, metric::Chebyshev<double>()), 7.625814583693912));
+    // BOOST_TEST(metric::entropy(v5, 3, 2.0, metric::Chebyshev<double>()) == 7.625814583693912);
+    // BOOST_TEST(metric::entropy(v6, 3, 2.0, metric::Chebyshev<double>()) == 7.625814583693912);
+    BOOST_TEST(float_eq(metric::entropy(v7, 3, 2.0, metric::Chebyshev<double>()), 7.625814583693912));
+
+    BOOST_TEST(float_eq(metric::entropy(v1),7.8197919635268081 ));
+    BOOST_TEST(float_eq(metric::entropy(v2), 7.8197919635268081));
+    BOOST_TEST(float_eq(metric::entropy(v3), 7.8197919635268081));
+    BOOST_TEST(float_eq(metric::entropy(v4), 7.8197919635268081));
+    BOOST_TEST(float_eq(metric::entropy(v7), 7.8197919635268081));
+
+    std::vector<std::string> v8 = { "AAA", "HJGJHFG", "BBB", "AAAA", "long long long long long long string" };
+    BOOST_TEST(float_eq(metric::entropy(v8, 3, 2.0, metric::Edit<int>()), 13.837977646415816));
+}
+
+BOOST_AUTO_TEST_CASE(mutual_information_float) {
+
+    std::vector<std::vector<double>> v11 = { { 5, 5 }, { 2, 2 }, { 3, 3 }, { 5, 5 } };
+    std::vector<std::vector<double>> v12 = { { 5, 5 }, { 2, 2 }, { 3, 3 }, { 1, 1 } };
+    metric::mutualInformation(v11, v12);
+
+    std::vector<std::deque<double>> v21 = { { 5, 5 }, { 2, 2 }, { 3, 3 }, { 5, 5 } };
+    std::vector<std::deque<double>> v22 = { { 5, 5 }, { 2, 2 }, { 3, 3 }, { 1, 1 } };
+    metric::mutualInformation(v21, v22);
+
+    std::deque<std::vector<double>> v31 = { { 5, 5 }, { 2, 2 }, { 3, 3 }, { 5, 5 } };
+    std::deque<std::vector<double>> v32 = { { 5, 5 }, { 2, 2 }, { 3, 3 }, { 1, 1 } };
+    metric::mutualInformation(v31, v32);
+
+    std::deque<std::deque<double>> v41 = { { 5, 5 }, { 2, 2 }, { 3, 3 }, { 5, 5 } };
+    std::deque<std::deque<double>> v42 = { { 5, 5 }, { 2, 2 }, { 3, 3 }, { 1, 1 } };
+    metric::mutualInformation(v41, v42);
+
+}
+BOOST_AUTO_TEST_CASE(mutual_information_int)
+{
+    /**
+     *  This test breaks compilation because of overload of the entropy method for integers is not implemented
+     *
+     */
+    // std::vector<std::vector<int>> v11 = { { 5, 5 }, { 2, 2 }, { 3, 3 }, { 5, 5 } };
+    // std::vector<std::vector<int>> v12 = { { 5, 5 }, { 2, 2 }, { 3, 3 }, { 1, 1 } };
+    // metric::mutualInformation(v11, v12);
+
+    // std::vector<std::deque<int>> v21 = { { 5, 5 }, { 2, 2 }, { 3, 3 }, { 5, 5 } };
+    // std::vector<std::deque<int>> v22 = { { 5, 5 }, { 2, 2 }, { 3, 3 }, { 1, 1 } };
+    // metric::mutualInformation(v21, v22);
+
+    // std::deque<std::vector<int>> v31 = { { 5, 5 }, { 2, 2 }, { 3, 3 }, { 5, 5 } };
+    // std::deque<std::vector<int>> v32 = { { 5, 5 }, { 2, 2 }, { 3, 3 }, { 1, 1 } };
+    // metric::mutualInformation(v31, v32);
+
+    // std::deque<std::deque<int>> v41 = { { 5, 5 }, { 2, 2 }, { 3, 3 }, { 5, 5 } };
+    // std::deque<std::deque<int>> v42 = { { 5, 5 }, { 2, 2 }, { 3, 3 }, { 1, 1 } };
+    // metric::mutualInformation(v41, v42);
+
+}
+BOOST_AUTO_TEST_CASE(variation_information) {
+    std::vector<std::vector<double>> v11 = { { 5, 5 }, { 2, 2 }, { 3, 3 }, { 5, 5 } };
+    std::vector<std::vector<double>> v12 = { { 5, 5 }, { 2, 2 }, { 3, 3 }, { 1, 1 } };
+    metric::variationOfInformation(v11, v12);
+    metric::variationOfInformation<decltype(v11), metric::Chebyshev<double>>(v11, v12);
+
+    std::vector<std::deque<double>> v21 = { { 5, 5 }, { 2, 2 }, { 3, 3 }, { 5, 5 } };
+    std::vector<std::deque<double>> v22 = { { 5, 5 }, { 2, 2 }, { 3, 3 }, { 1, 1 } };
+    metric::variationOfInformation(v21, v22);
+    metric::variationOfInformation<decltype(v21), metric::Chebyshev<double>>(v21, v22);
+
+    std::deque<std::vector<double>> v31 = { { 5, 5 }, { 2, 2 }, { 3, 3 }, { 5, 5 } };
+    std::deque<std::vector<double>> v32 = { { 5, 5 }, { 2, 2 }, { 3, 3 }, { 1, 1 } };
+    metric::variationOfInformation(v31, v32);
+    metric::variationOfInformation<decltype(v31), metric::Chebyshev<double>>(v31, v32);
+
+    std::deque<std::deque<double>> v41 = { { 5, 5 }, { 2, 2 }, { 3, 3 }, { 5, 5 } };
+    std::deque<std::deque<double>> v42 = { { 5, 5 }, { 2, 2 }, { 3, 3 }, { 1, 1 } };
+    metric::variationOfInformation(v41, v42);
+    metric::variationOfInformation<decltype(v41), metric::Chebyshev<double>>(v41, v42);
+}
+
+BOOST_AUTO_TEST_CASE(variation_information_normalized)
+{
+    std::vector<std::vector<double>> v11 = { { 5, 5 }, { 2, 2 }, { 3, 3 }, { 5, 5 } };
+    std::vector<std::vector<double>> v12 = { { 5, 5 }, { 2, 2 }, { 3, 3 }, { 1, 1 } };
+    metric::variationOfInformation_normalized(v11, v12);
+
+    std::vector<std::deque<double>> v21 = { { 5, 5 }, { 2, 2 }, { 3, 3 }, { 5, 5 } };
+    std::vector<std::deque<double>> v22 = { { 5, 5 }, { 2, 2 }, { 3, 3 }, { 1, 1 } };
+    metric::variationOfInformation_normalized(v21, v22);
+
+    std::deque<std::vector<double>> v31 = { { 5, 5 }, { 2, 2 }, { 3, 3 }, { 5, 5 } };
+    std::deque<std::vector<double>> v32 = { { 5, 5 }, { 2, 2 }, { 3, 3 }, { 1, 1 } };
+    metric::variationOfInformation_normalized(v31, v32);
+
+    std::deque<std::deque<double>> v41 = { { 5, 5 }, { 2, 2 }, { 3, 3 }, { 5, 5 } };
+    std::deque<std::deque<double>> v42 = { { 5, 5 }, { 2, 2 }, { 3, 3 }, { 1, 1 } };
+    metric::variationOfInformation_normalized(v41, v42);
+
+}
+
+BOOST_AUTO_TEST_CASE(VOI) {
+
+    std::vector<std::vector<double>> v11 = { { 5, 5 }, { 2, 2 }, { 3, 3 }, { 5, 5 } };
+    std::vector<std::vector<double>> v12 = { { 5, 5 }, { 2, 2 }, { 3, 3 }, { 1, 1 } };
+    metric::VOI<long double>()(v11, v12);
+    metric::VOI_normalized<long double>()(v11, v12);
+    metric::VOI_kl<long double>()(v11, v12);
+    metric::VOI_normalized_kl<long double>()(v11, v12);
+    
+    std::vector<std::deque<double>> v21 = { { 5, 5 }, { 2, 2 }, { 3, 3 }, { 5, 5 } };
+    std::vector<std::deque<double>> v22 = { { 5, 5 }, { 2, 2 }, { 3, 3 }, { 1, 1 } };
+    metric::VOI<long double>()(v21, v22);
+    metric::VOI_normalized<long double>()(v21, v22);
+    metric::VOI_kl<long double>()(v21, v22);
+    metric::VOI_normalized_kl<long double>()(v21, v22);
+
+    std::deque<std::vector<double>> v31 = { { 5, 5 }, { 2, 2 }, { 3, 3 }, { 5, 5 } };
+    std::deque<std::vector<double>> v32 = { { 5, 5 }, { 2, 2 }, { 3, 3 }, { 1, 1 } };
+    metric::VOI<long double>()(v31, v32);
+    metric::VOI_normalized<long double>()(v31, v32);
+    metric::VOI_kl<long double>()(v31, v32);
+    metric::VOI_normalized_kl<long double>()(v31, v32);
+
+    std::deque<std::deque<double>> v41 = { { 5, 5 }, { 2, 2 }, { 3, 3 }, { 5, 5 } };
+    std::deque<std::deque<double>> v42 = { { 5, 5 }, { 2, 2 }, { 3, 3 }, { 1, 1 } };
+    metric::VOI<long double>()(v41, v42);
+    metric::VOI_normalized<long double>()(v41, v42);
+    metric::VOI_kl<long double>()(v41, v42);
+    metric::VOI_normalized_kl<long double>()(v41, v42);
+}

From 56f76006867d927493c1496a2b4f2b95434a90cd Mon Sep 17 00:00:00 2001
From: LeilaPiri <54840627+LeilaPiri@users.noreply.github.com>
Date: Fri, 20 Mar 2020 23:13:51 +0100
Subject: [PATCH 2/4] Reverse diffusion algorithm (#125)

* adding to repository on behalf of another programmer

* Update modules/utils/graph/Redif.cpp

* Apply suggestions from code review

* apply suggestions from CR, require testing

+ fix missing const&
+ remove unneeded template parameters
+ simplify passing metric by value

Co-authored-by: Pauls <pauls@frigate.frigate.lv>
Co-authored-by: jura-gresko <57671564+jura-gresko@users.noreply.github.com>
Co-authored-by: jura-gresko <juragresko@gmail.com>
---
 modules/mapping/Redif.cpp | 326 ++++++++++++++++++++++++++++++++++++++
 modules/mapping/Redif.hpp |  60 +++++++
 2 files changed, 386 insertions(+)
 create mode 100644 modules/mapping/Redif.cpp
 create mode 100644 modules/mapping/Redif.hpp

diff --git a/modules/mapping/Redif.cpp b/modules/mapping/Redif.cpp
new file mode 100644
index 00000000..33adf471
--- /dev/null
+++ b/modules/mapping/Redif.cpp
@@ -0,0 +1,326 @@
+#include "Redif.hpp"
+
+
+namespace metric
+{
+	template <typename Tv>
+	template <class Metric>
+	Redif<Tv>::Redif(
+	    const std::vector<std::vector<Tv>>& trainData,
+	    size_t nNeighbors /* = 10 */,
+	    size_t nIter /* = 15 */,
+	    Metric metric = Metric()
+    )
+	{
+		this->nNeighbors = nNeighbors;
+		xTrain = blaze::DynamicMatrix<Tv>(trainData.size(), trainData[0].size());
+		for (size_t i = 0; i < xTrain.rows(); i++)
+			for (size_t j = 0; j < xTrain.columns(); j++)
+				xTrain(i, j) = trainData[i][j];
+		trainModel(nIter, metric);
+	}
+
+
+	template <typename Tv>
+	template <class Metric>
+	void Redif<Tv>::trainModel(size_t nIter, Metric metric)
+	{
+		blaze::DynamicMatrix<Tv> xTrainDe(this->xTrain);
+		const size_t n = xTrainDe.rows();
+
+		/// Iterative backward Diffusion process
+		for (size_t i = 0; i < nIter; i++)
+		{
+			blaze::DynamicMatrix<Tv> localDistMat = getLocalDistMatrix(xTrainDe, metric);
+			blaze::DynamicMatrix<Tv> graphLap = calcWeightedGraphLaplacian(localDistMat);
+			this->LArray.push_back(graphLap);
+			blaze::DynamicMatrix<Tv> D(n, n);
+			for (size_t i = 0; i < D.rows(); i++)
+				D(i, i) = graphLap(i, i);
+			blaze::DynamicMatrix<Tv> invLeft(n, n);
+			invLeft = blaze::inv(D + 0.25 * graphLap);
+			xTrainDe = invLeft * (D * xTrainDe);
+		}
+
+		return xTrainDe;
+	}
+
+	template <typename Tv>
+	template <class Metric>
+	blaze::DynamicMatrix<Tv> Redif<Tv>::getLocalDistMatrix(
+	    const blaze::DynamicMatrix<Tv>& dataSample,
+	    Metric metric
+	)
+	{
+		size_t n = dataSample.rows();
+
+		blaze::DynamicMatrix<size_t> ixx(nNeighbors, n);
+		for (size_t i = 0; i < nNeighbors; i++)
+			for (size_t j = 0; j < n; j++)
+				ixx(i, j) = j;
+
+		blaze::DynamicMatrix<Tv> dist(n, n);
+		for (size_t i = 0; i < n; i++)
+			for (size_t j = 0; j < n; j++)
+			{
+				if (i == j)
+				{
+					dist(i, j) = 0;
+					continue;
+				}
+				blaze::DynamicVector<Tv, rowVector> rowi = row(dataSample, i);
+				blaze::DynamicVector<Tv, rowVector> rowj = row(dataSample, j);
+
+
+				std::vector<Tv> veci, vecj;
+				for (size_t k = 0; k < rowi.size(); ++k)
+				{
+					veci.push_back(rowi[k]);
+					vecj.push_back(rowj[k]);
+				}				
+
+				dist(i, j) = metric(veci, vecj);
+			}
+		blaze::DynamicMatrix<size_t> knnMat(nNeighors, n);
+		blaze::DynamicMatrix<Tv> knnDistMat(nNeighbors, n);
+
+		for (size_t i = 0; i < n; i++)
+		{
+			blaze::DynamicVector<size_t> tempIdx(n);
+			blaze::DynamicVector<Tv> tempDist(n);
+			for (size_t j = 0; j < n; j++)
+			{
+				tempIdx[j] = j;
+				tempDist[j] = dist(j, i);
+			}
+			Quicksort(tempIdx, tempDist, 0, n - 1);
+			for (size_t j = 0; j < nNeighbors; j++)
+			{
+				knnMat(j, i) = tempIdx[j + 1];
+				knnDistMat(j, i) = tempDist[j + 1];
+			}
+		}
+
+		blaze::DynamicMatrix<Tv> localDist(n, n);
+
+		for (size_t i = 0; i < nNeighbors; i++)
+			for (size_t j = 0; j < n; j++)
+				localDist(ixx(i, j), knnMat(i, j)) = knnDistMat(i, j);
+
+		blaze::DynamicMatrix<Tv> retLocalDist(n, n);
+		for (size_t i = 0; i < n; i++)
+			for (size_t j = 0; j < n; j++)
+				retLocalDist(i, j) = (localDist(i, j) > localDist(j, i)) ? localDist(i, j) : localDist(j, i);
+
+		return retLocalDist;
+	}
+
+	template <typename Tv>
+	blaze::DynamicMatrix<Tv> Redif<Tv>::calcWeightedGraphLaplacian(const blaze::DynamicMatrix<Tv>& localDist)
+
+	{
+		size_t n = localDist.rows();
+
+		/// Degree matrix
+		blaze::DynamicMatrix<Tv> D(n, n);
+		for (size_t i = 0; i < D.rows(); i++)
+		{
+			double sum = 0;
+			for (size_t j = 0; j < n; j++)
+				sum += localDist(i, j);
+			if (sum == 0) sum = 1;
+			D(i, i) = 1 / sum;
+		}
+
+		blaze::DynamicMatrix<Tv> A(localDist.rows(), localDist.columns());
+		A = D * localDist * D / n;
+
+		blaze::DynamicMatrix<Tv> D1(n, n);
+		for (size_t i = 0; i < D.rows(); i++)
+		{
+			double sum = 0;
+			for (size_t j = 0; j < n; j++)
+				sum += A(i, j);
+			if (sum == 0) sum = 1 / n;
+			D1(i, i) = sum;
+		}
+
+		blaze::DynamicMatrix<Tv> graphLap(n, n);
+		graphLap = D1 - A;
+
+		return graphLap;
+	}
+
+	template <typename Tv>
+	template <class Metric>
+	blaze::DynamicMatrix<Tv> Redif<Tv>::encode(
+	    const blaze::DynamicMatrix<Tv>& x,
+	    blaze::DynamicVector<Tv>& l_idx,
+	    Metric metric = Metric()
+    )
+	{
+		size_t nTrain = xTrain.rows();
+		size_t nX = x.rows();
+		size_t nIter = LArray.size();
+		/// Find nearest Neighbors of each record of x in xTrain
+		size_t l = 1; /// Number of nearest neighbors
+		l_idx = blaze::DynamicVector<Tv>(nX);
+		unsigned size_t k_check = (unsigned size_t)sqrt(nTrain); /// Parameter for checking model on the record
+
+		/// Compute distances
+		for (size_t i = 0; i < nX; i++)
+		{
+			blaze::DynamicMatrix<Tv> distTotal(nTrain + 1, nTrain + 1);
+			for (size_t k = 0; k < nTrain + 1; k++)
+				for (size_t j = 0; j < nTrain + 1; j++)
+				{
+					if (k == j) continue;
+
+					blaze::DynamicVector<Tv, rowVector> rowk;
+					blaze::DynamicVector<Tv, rowVector> rowj;
+
+					if (k < nTrain) rowk = row(xTrain, k);
+					else rowk = row(x, i);
+					if (j < nTrain) rowj = row(xTrain, j);
+					else rowj = row(x, i);
+
+					std::vector<Tv> veck, vecj;
+					for (size_t s = 0; s < rowk.size(); ++s)
+					{
+						veck.push_back(rowk[s]);
+						vecj.push_back(rowj[s]);
+					}
+
+					distTotal(k, j) = metric(veck, vecj);
+				}
+
+			blaze::DynamicVector<Tv, columnVector> dist = column(distTotal, distTotal.columns() - 1);
+
+			blaze::DynamicVector<size_t> idxVec(dist.size());
+			for (size_t k = 0; k < idxVec.size(); k++) idxVec[k] = k;
+			Quicksort(idxVec, dist, 0, idxVec.size() - 1);
+
+			double oneDist = 0;
+			if (dist[0] == dist[1])
+			{
+				size_t j = 0;
+				j = (idxVec[0] > idxVec[1]) ? 0 : 1;
+				if (j == 0)
+				{
+					l_idx[i] = idxVec[1];
+					oneDist = dist[1];
+				}
+				else
+				{
+					l_idx[i] = idxVec[0];
+					oneDist = dist[0];
+				}
+			}
+			else
+			{
+				l_idx[i] = idxVec[1];
+				oneDist = dist[1];
+			}
+			/// Check if Model is able to denoise Test Posize_t
+			blaze::DynamicVector<Tv, columnVector> temp = column(distTotal, l_idx[i]);
+			for (size_t k = 0; k < idxVec.size(); k++) idxVec[k] = k;
+			Quicksort(idxVec, temp, 0, idxVec.size() - 1);
+			if (oneDist > temp[k_check - 1])
+				cout << "Error. Testposize_t of Index " << i << " cannot be denoised based on underlying model.\n";
+		}
+
+		blaze::DynamicMatrix<Tv> xEncoded = xTrain;
+		/// Iterative backward Diffusion process
+		for (size_t i = 0; i < nX; i++)
+			for (size_t j = 0; j < xEncoded.columns(); j++)
+				xEncoded(l_idx[i], j) = x(i, j);
+
+		for (size_t i = 0; i < nIter; i++)
+		{
+			blaze::DynamicMatrix<Tv> D(nTrain, nTrain);
+			for (size_t k = 0; k < D.rows(); k++)
+				D(k, k) = LArray[i](k, k);
+			xEncoded = blaze::inv(D + 0.25*LArray[i])*(D*xEncoded);
+		}
+
+		return xEncoded;
+	}
+
+	template <typename Tv>
+	blaze::DynamicMatrix<Tv> Redif<Tv>::decode(const blaze::DynamicMatrix<Tv>& xEncoded, const blaze::DynamicMatrix<Tv>& xTrainEncoded, const blaze::DynamicVector<size_t>& l_idx)
+
+	{
+		size_t n = xTrainEncoded.rows();
+		size_t nIter = LArray.size();
+		for (size_t i = 0; i < l_idx.size(); i++)
+			for (size_t j = 0; j < xTrainEncoded.columns(); j++)
+				xTrainEncoded(l_idx[i], j) = xEncoded(i, j);
+
+		for (size_t i = 0; i < nIter; i++)
+		{
+			blaze::DynamicMatrix<Tv> temp = LArray[nIter - i - 1];
+			blaze::DynamicMatrix<Tv> D(n, n);
+			for (size_t k = 0; k < n; k++)
+				D(k, k) = temp(k, k);
+			xTrainEncoded = blaze::inv(D)*(D + 0.25*LArray[nIter - i - 1])*xTrainEncoded;
+		}
+		blaze::DynamicMatrix<Tv> xDecoded = blaze::DynamicMatrix<Tv>(l_idx.size(), xTrainEncoded.columns());
+		for (size_t i = 0; i < xDecoded.rows(); i++)
+			for (size_t j = 0; j < xDecoded.columns(); j++)
+				xDecoded(i, j) = xTrainEncoded(l_idx[i], j);
+
+		return xDecoded;
+	}
+
+	template <typename Tv>
+	template<bool flag>
+	size_t Redif<Tv>::partition(blaze::DynamicVector<size_t> & idx, blaze::DynamicVector<Tv, flag> & val, size_t start, size_t end)
+	{
+		double pivot = val[end];
+
+		size_t P_index = start;
+		size_t i;
+		size_t idxTemp;
+		double valTemp;
+
+		for (i = start; i < end; i++)
+		{
+			if (val[i] <= pivot)
+			{
+				valTemp = val[i];
+				val[i] = val[P_index];
+				val[P_index] = valTemp;
+
+				idxTemp = idx[i];
+				idx[i] = idx[P_index];
+				idx[P_index] = idxTemp;
+
+				P_index++;
+			}
+		}
+
+		valTemp = val[end];
+		val[end] = val[P_index];
+		val[P_index] = valTemp;
+
+		idxTemp = idx[end];
+		idx[end] = idx[P_index];
+		idx[P_index] = idxTemp;
+
+		return P_index;
+	}
+
+	template <typename Tv>
+	template<bool flag>
+	void Redif<Metric, Tv>::Quicksort(blaze::DynamicVector<size_t> & idx, blaze::DynamicVector<Tv, flag> & data, size_t start, size_t end)
+	{
+		if (start < end)
+		{
+			int P_index = partition(idx, data, start, end);
+			if (P_index != 0) Quicksort(idx, data, start, P_index - 1);
+			Quicksort(idx, data, P_index + 1, end);
+		}
+	}
+
+
+}
diff --git a/modules/mapping/Redif.hpp b/modules/mapping/Redif.hpp
new file mode 100644
index 00000000..979a6d1f
--- /dev/null
+++ b/modules/mapping/Redif.hpp
@@ -0,0 +1,60 @@
+#pragma once
+
+#include "3rdparty/blaze/Blaze.h"
+#include "modules/distance.hpp"
+
+#include <vector>
+
+namespace metric
+{
+	template<typename Tv = double>
+	class Redif
+	{
+	public:
+	    template <class Metric = metric::Euclidian<Tv>>
+		Redif(
+		    const std::vector<std::vector<Tv>>& trainData,
+		    size_t nNeighbors = 10,
+		    size_t nIter = 15
+		    Metric metric = Metric()
+        );
+		~Redif() = default;
+
+		template <class Metric = metric::Euclidian<Tv>>
+		blaze::DynamicMatrix<Tv> encode(
+		    const blaze::DynamicMatrix<Tv>& x,
+		    blaze::DynamicVector<Tv> & l_idx,
+		    Metric metric = Metric()
+        );
+		blaze::DynamicMatrix<Tv> decode(
+		    const blaze::DynamicMatrix<Tv>& xEncoded,
+		    const blaze::DynamicMatrix<Tv>& xTrainEncoded,
+		    const blaze::DynamicVector<size_t>& l_idx
+        );
+	private:
+
+		template<bool flag = blaze::rowVector>
+		void Quicksort(
+		    blaze::DynamicVector<size_t> & idx,
+		    blaze::DynamicVector<Tv, flag> & data,
+		    size_t start,
+		    size_t end
+        );
+
+		template<bool flag = blaze::rowVector>
+		size_t partition(
+		    blaze::DynamicVector<size_t> & idx,
+		    blaze::DynamicVector<Tv, flag> & data,
+		    size_t start,
+		    size_t end
+        );
+
+		blaze::DynamicMatrix<Tv> getLocalDistMatrix(const blaze::DynamicMatrix<Tv>& dataSample);
+		blaze::DynamicMatrix<Tv> calcWeightedGraphLaplacian(const blaze::DynamicMatrix<Tv>& localDist);
+		blaze::DynamicMatrix<Tv> trainModel(size_t nIter);
+
+		blaze::DynamicMatrix<Tv> xTrain;
+		std::vector<blaze::DynamicMatrix<Tv>> LArray;
+		size_t nNeighbors;
+	};
+}

From 05838129fc615f9a6a455441a31146b86073dd76 Mon Sep 17 00:00:00 2001
From: Stepka <micromillioner@yahoo.com>
Date: Sat, 21 Mar 2020 00:21:29 +0200
Subject: [PATCH 3/4] K random metrics (#107)

* [Kolmogorov-Smirnov] added implementation for Kolmogorov-Smirnov distance

* [Kolmogorov-Smirnov] removed unneccesary files

* [Kolmogorov-Smirnov] Updated Discrete distribution to be created from samples, updated mean, variance and cdf.

* little fix

* [Kolmogorov-Smirnov distance] Updated distance to use pmq's cdf

* [Random EMD] added earth's mover distance based on cdf

* Added Cramer-von Mises distance

* Added tests for some k-random distances

* Added some documentation

* [K-random_metrics] Updated PMQ constructor - create distributions from data histogram

Co-authored-by: Michael Welsch <welsch@panda.technology>
---
 examples/distance_examples/CMakeLists.txt     |  10 +-
 .../cramervon_mises_distance_example.cpp      |  62 +++++++++
 .../kolmogorov_smirnov_distance_example.cpp   |  40 ++++++
 .../distance_examples/random_emd_example.cpp  |  76 +++++++++++
 modules/distance.hpp                          |   4 +
 modules/distance/k-random/CramervonNises.cpp  |  50 ++++++++
 modules/distance/k-random/CramervonNises.hpp  |  63 +++++++++
 .../distance/k-random/KolmogorovSmirnov.cpp   |  50 ++++++++
 .../distance/k-random/KolmogorovSmirnov.hpp   |  50 ++++++++
 modules/distance/k-random/RandomEMD.cpp       |  51 ++++++++
 modules/distance/k-random/RandomEMD.hpp       |  64 +++++++++
 modules/mapping/KOC.hpp                       | 121 +++++++++---------
 modules/utils/poor_mans_quantum.hpp           |  10 +-
 .../distributions/Discrete.cpp                |  36 ++++--
 .../poor_mans_quantum/poor_mans_quantum.cpp   |  65 ++++++----
 tests/distance_tests/CMakeLists.txt           |  72 ++++++++---
 tests/distance_tests/README.md                |  30 +++++
 .../distance_tests/cramervon_mises_tests.cpp  |  92 +++++++++++++
 .../kolmogorov_smirnov_tests.cpp              |  95 ++++++++++++++
 tests/distance_tests/random_emd_tests.cpp     |  91 +++++++++++++
 20 files changed, 1011 insertions(+), 121 deletions(-)
 create mode 100644 examples/distance_examples/cramervon_mises_distance_example.cpp
 create mode 100644 examples/distance_examples/kolmogorov_smirnov_distance_example.cpp
 create mode 100644 examples/distance_examples/random_emd_example.cpp
 create mode 100644 modules/distance/k-random/CramervonNises.cpp
 create mode 100644 modules/distance/k-random/CramervonNises.hpp
 create mode 100644 modules/distance/k-random/KolmogorovSmirnov.cpp
 create mode 100644 modules/distance/k-random/KolmogorovSmirnov.hpp
 create mode 100644 modules/distance/k-random/RandomEMD.cpp
 create mode 100644 modules/distance/k-random/RandomEMD.hpp
 create mode 100644 tests/distance_tests/README.md
 create mode 100644 tests/distance_tests/cramervon_mises_tests.cpp
 create mode 100644 tests/distance_tests/kolmogorov_smirnov_tests.cpp
 create mode 100644 tests/distance_tests/random_emd_tests.cpp

diff --git a/examples/distance_examples/CMakeLists.txt b/examples/distance_examples/CMakeLists.txt
index 59115998..0350f1be 100644
--- a/examples/distance_examples/CMakeLists.txt
+++ b/examples/distance_examples/CMakeLists.txt
@@ -5,18 +5,18 @@ project(distance_examples)
 set(Boost_DEBUG on)  
 find_package(Boost)
 
-if(UNIX)
+if(UNIX) 
     find_package(Threads REQUIRED)
 endif(UNIX)
 
 include_directories(
     ${PROJECT_SOURCE_DIR}
     ${PROJECT_SOURCE_DIR}/../../3rdparty/jpeglib 
-    )
-
-include_directories(${Boost_INCLUDE_DIRS}) 
+    ) 
 	 
-file(COPY "assets" DESTINATION ".")
+include_directories(${Boost_INCLUDE_DIRS})  
+	
+file(COPY "assets" DESTINATION ".") 
 
 file(GLOB EXAMPLE_SRCS RELATIVE ${CMAKE_CURRENT_SOURCE_DIR} *.cpp) 
 find_package(JPEG) 
diff --git a/examples/distance_examples/cramervon_mises_distance_example.cpp b/examples/distance_examples/cramervon_mises_distance_example.cpp
new file mode 100644
index 00000000..3ecf3bac
--- /dev/null
+++ b/examples/distance_examples/cramervon_mises_distance_example.cpp
@@ -0,0 +1,62 @@
+/*
+This Source Code Form is subject to the terms of the Mozilla Public
+License, v. 2.0. If a copy of the MPL was not distributed with this
+file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+Copyright (c) 2019 Panda Team
+*/
+
+#include <vector>
+#include <iostream>
+#include <chrono>
+#include "../../modules/distance.hpp"
+#include "assets/test_data.cpp"
+
+
+int main()
+{
+	    /******************** examples for Cramer-von Nises Distance **************************/
+    // example for picture
+    std::cout << "Cramer-von Nises distance example have started" << std::endl;
+    std::cout << "" << std::endl;	
+
+	/*** here are some data records ***/
+	std::vector<double> samples_1 = { 0, 1, 2, 3, 3, 2, 1, 0, 2, 2 };
+	std::vector<double> samples_2 = { 0, 0, 2, 3, 3, 2, 1, 0, 2, 2 };
+	   
+    metric::CramervonNises<std::vector<double>, double> distance_1;
+
+    auto t1 = std::chrono::steady_clock::now();
+    auto result = distance_1(samples_1, samples_2);
+    auto t2 = std::chrono::steady_clock::now();
+    std::cout << "result: " << result
+              << " (Time = " << double(std::chrono::duration_cast<std::chrono::microseconds>(t2 - t1).count()) / 1000000
+              << " s)" << std::endl;
+    std::cout << "" << std::endl;
+
+	//
+	   
+    metric::CramervonNises<std::vector<double>, double> distance_2(0.0001);
+
+    t1 = std::chrono::steady_clock::now();
+    result = distance_2(samples_1, samples_2);
+    t2 = std::chrono::steady_clock::now();
+    std::cout << "result: " << result
+              << " (Time = " << double(std::chrono::duration_cast<std::chrono::microseconds>(t2 - t1).count()) / 1000000
+              << " s)" << std::endl;
+    std::cout << "" << std::endl;
+
+	//
+	   
+    metric::CramervonNises<std::vector<double>, double> distance_3(0.1);
+
+    t1 = std::chrono::steady_clock::now();
+    result = distance_3(samples_1, samples_2);
+    t2 = std::chrono::steady_clock::now();
+    std::cout << "result: " << result
+              << " (Time = " << double(std::chrono::duration_cast<std::chrono::microseconds>(t2 - t1).count()) / 1000000
+              << " s)" << std::endl;
+    std::cout << "" << std::endl;
+
+    return 0;
+}
diff --git a/examples/distance_examples/kolmogorov_smirnov_distance_example.cpp b/examples/distance_examples/kolmogorov_smirnov_distance_example.cpp
new file mode 100644
index 00000000..46be3104
--- /dev/null
+++ b/examples/distance_examples/kolmogorov_smirnov_distance_example.cpp
@@ -0,0 +1,40 @@
+/*
+This Source Code Form is subject to the terms of the Mozilla Public
+License, v. 2.0. If a copy of the MPL was not distributed with this
+file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+Copyright (c) 2019 Panda Team
+*/
+
+#include <vector>
+#include <iostream>
+#include <chrono>
+#include "../../modules/distance.hpp"
+#include "assets/test_data.cpp"
+
+
+int main()
+{
+	    /******************** examples for Kolmogorov-Smirnov Distance **************************/
+    // example for picture
+    std::cout << "Kolmogorov-Smirnov distance example have started" << std::endl;
+    std::cout << "" << std::endl;	
+
+	/*** here are some data records ***/
+	std::vector<double> samples_1 = { 0, 1, 2, 3, 3, 2, 1, 0, 2, 2 };
+	std::vector<double> samples_2 = { 0, 0, 2, 3, 3, 2, 1, 0, 2, 2 };
+	//std::vector<double> samples_1 = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 };
+	//std::vector<double> samples_2 = { 1, 2, 3, 4, 5 };
+	//   
+    metric::KolmogorovSmirnov<std::vector<double>, double> distance;
+
+    auto t1 = std::chrono::steady_clock::now();
+    auto result1 = distance(samples_1, samples_2);
+    auto t2 = std::chrono::steady_clock::now();
+    std::cout << "result: " << result1
+              << " (Time = " << double(std::chrono::duration_cast<std::chrono::microseconds>(t2 - t1).count()) / 1000000
+              << " s)" << std::endl;
+    std::cout << "" << std::endl;
+
+    return 0;
+}
diff --git a/examples/distance_examples/random_emd_example.cpp b/examples/distance_examples/random_emd_example.cpp
new file mode 100644
index 00000000..df95ac92
--- /dev/null
+++ b/examples/distance_examples/random_emd_example.cpp
@@ -0,0 +1,76 @@
+/*
+This Source Code Form is subject to the terms of the Mozilla Public
+License, v. 2.0. If a copy of the MPL was not distributed with this
+file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+Copyright (c) 2019 Panda Team
+*/
+
+#include <vector>
+#include <iostream>
+#include <chrono>
+#include "../../modules/distance.hpp"
+#include "assets/test_data.cpp"
+
+
+int main()
+{
+	    /******************** examples for Earth Mover’s Distance **************************/
+    // example for picture
+    std::cout << "Earth Mover's distance example have started" << std::endl;
+    std::cout << "" << std::endl;	
+
+	/*** here are some data records ***/
+	std::vector<double> samples_1 = { 0, 1, 2, 3, 3, 2, 1, 0, 2, 2 };
+	std::vector<double> samples_2 = { 0, 0, 2, 3, 3, 2, 1, 0, 2, 2 };
+	   
+    metric::RandomEMD<std::vector<double>, double> distance_1;
+
+
+	//
+ //   typedef int emd_Type;
+
+ //   auto cost_mat = metric::EMD_details::ground_distance_matrix_of_2dgrid<emd_Type>(3, 3);
+ //   auto maxCost = metric::EMD_details::max_in_distance_matrix(cost_mat);
+
+ //   metric::EMD<emd_Type> distance_orig(cost_mat, maxCost);
+
+	//print_matrix(cost_mat);
+
+
+
+
+    auto t1 = std::chrono::steady_clock::now();
+    auto result = distance_1(samples_1, samples_2);
+    auto t2 = std::chrono::steady_clock::now();
+    std::cout << "result: " << result
+              << " (Time = " << double(std::chrono::duration_cast<std::chrono::microseconds>(t2 - t1).count()) / 1000000
+              << " s)" << std::endl;
+    std::cout << "" << std::endl;
+
+	//
+	   
+    metric::RandomEMD<std::vector<double>, double> distance_2(0.0001);
+
+    t1 = std::chrono::steady_clock::now();
+    result = distance_2(samples_1, samples_2);
+    t2 = std::chrono::steady_clock::now();
+    std::cout << "result: " << result
+              << " (Time = " << double(std::chrono::duration_cast<std::chrono::microseconds>(t2 - t1).count()) / 1000000
+              << " s)" << std::endl;
+    std::cout << "" << std::endl;
+
+	//
+	   
+    metric::RandomEMD<std::vector<double>, double> distance_3(0.1);
+
+    t1 = std::chrono::steady_clock::now();
+    result = distance_3(samples_1, samples_2);
+    t2 = std::chrono::steady_clock::now();
+    std::cout << "result: " << result
+              << " (Time = " << double(std::chrono::duration_cast<std::chrono::microseconds>(t2 - t1).count()) / 1000000
+              << " s)" << std::endl;
+    std::cout << "" << std::endl;
+
+    return 0;
+}
diff --git a/modules/distance.hpp b/modules/distance.hpp
index 59f7305a..afd1ea21 100644
--- a/modules/distance.hpp
+++ b/modules/distance.hpp
@@ -32,6 +32,10 @@ Copyright (c) PANDA Team
 #include "distance/k-structured/kohonen_distance.hpp"
 
 #include "distance/k-related/L1.hpp"
+
 #include "distance/k-random/VOI.hpp"
+#include "distance/k-random/KolmogorovSmirnov.hpp"
+#include "distance/k-random/RandomEMD.hpp"
+#include "distance/k-random/CramervonNises.hpp"
 
 #endif //_METRIC_DISTANCE_HPP
diff --git a/modules/distance/k-random/CramervonNises.cpp b/modules/distance/k-random/CramervonNises.cpp
new file mode 100644
index 00000000..30e269af
--- /dev/null
+++ b/modules/distance/k-random/CramervonNises.cpp
@@ -0,0 +1,50 @@
+/*
+This Source Code Form is subject to the terms of the Mozilla Public
+License, v. 2.0. If a copy of the MPL was not distributed with this
+file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+Copyright (c) 2020 Panda Team
+*/
+#ifndef _METRIC_DISTANCE_K_RANDOM_CRAMER_VON_NISES_CPP
+#define _METRIC_DISTANCE_K_RANDOM_CRAMER_VON_NISES_CPP
+
+#include "CramervonNises.hpp"
+#include "../../utils/poor_mans_quantum.hpp"
+
+namespace metric {
+
+template <typename Sample, typename D>
+auto CramervonNises<typename Sample, typename D>::operator()(const Sample& sample_1, const Sample& sample_2) -> distance_return_type
+{
+	PMQ pmq_1(sample_1);
+	PMQ pmq_2(sample_2);
+		
+	// find the sum of squared differences between the two cdfs
+
+	Sample concat_data;
+	for (int i = 0; i < sample_1.size(); i++)
+	{
+		concat_data.push_back(sample_1[i]);
+	}
+	for (int i = 0; i < sample_2.size(); i++)
+	{
+		concat_data.push_back(sample_2[i]);
+	}
+	std::sort(concat_data.begin(), concat_data.end());
+
+	//
+	
+	D area = 0;
+	double step = (concat_data[concat_data.size() - 1] - concat_data[0]) * precision_;
+
+	for (double value = concat_data[0]; value <= concat_data[concat_data.size() - 1]; value += step)
+	{
+		area += (pmq_1.cdf(value) - pmq_2.cdf(value)) * (pmq_1.cdf(value) - pmq_2.cdf(value)) * step;
+	}
+
+    return sqrt(area);
+}
+
+
+}  // namespace metric
+#endif
diff --git a/modules/distance/k-random/CramervonNises.hpp b/modules/distance/k-random/CramervonNises.hpp
new file mode 100644
index 00000000..c343afc3
--- /dev/null
+++ b/modules/distance/k-random/CramervonNises.hpp
@@ -0,0 +1,63 @@
+/*
+This Source Code Form is subject to the terms of the Mozilla Public
+License, v. 2.0. If a copy of the MPL was not distributed with this
+file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+Copyright (c) 2020 Panda Team
+*/
+
+#ifndef _METRIC_DISTANCE_K_RANDOM_CRAMER_VON_NISES_HPP
+#define _METRIC_DISTANCE_K_RANDOM_CRAMER_VON_NISES_HPP
+
+namespace metric {
+	
+
+/**
+ * @brief
+ *
+ * The Cramér-von Mises (CM) distance is obtained by summing the squared difference between the two empirical CDFs along the x-axis (and then taking the square root of the sum to make it an actual distance). 
+ * The relationship between the CM distance and the EMD is analogous to that of the L1 and L2 norms. 
+ * CDFs is calculated by PMQ.
+ * 
+ * @tparam Sample - sample type
+ * @tparam D - distance return type
+ */
+template <typename Sample, typename D = double>
+struct CramervonNises {
+    using distance_return_type = D;
+
+    /**
+     * @brief Construct a new Cramer-von Nises object
+     *
+     */
+    explicit CramervonNises()
+    {
+    }
+
+    /**
+     * @brief Construct a new Cramer-von Nises object
+     *
+     * @param precision  for integration. Should be in (0, 1). Less means more accurate.
+     */
+    explicit CramervonNises(double precision) : precision_(precision)
+    {
+    }
+
+    /**
+     * @brief calculate Cramer-von Nises distance between two samples
+     *
+     * @param sample_1 first sample
+     * @param sample_2 second sample
+     * @return distance
+     */
+    distance_return_type operator()(const Sample& sample_1, const Sample& sample_2);
+
+private:
+	double precision_ = 0.01;
+};
+
+
+}  // namespace metric
+
+#include "CramervonNises.cpp"
+#endif
diff --git a/modules/distance/k-random/KolmogorovSmirnov.cpp b/modules/distance/k-random/KolmogorovSmirnov.cpp
new file mode 100644
index 00000000..fe3931a5
--- /dev/null
+++ b/modules/distance/k-random/KolmogorovSmirnov.cpp
@@ -0,0 +1,50 @@
+/*
+This Source Code Form is subject to the terms of the Mozilla Public
+License, v. 2.0. If a copy of the MPL was not distributed with this
+file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+Copyright (c) 2020 Panda Team
+*/
+#ifndef _METRIC_DISTANCE_K_RANDOM_KOLMOGOROV_SMIRNOV_CPP
+#define _METRIC_DISTANCE_K_RANDOM_KOLMOGOROV_SMIRNOV_CPP
+
+#include "KolmogorovSmirnov.hpp"
+#include "../../utils/poor_mans_quantum.hpp"
+
+namespace metric {
+
+template <typename Sample, typename D>
+auto KolmogorovSmirnov<typename Sample, typename D>::operator()(const Sample& sample_1, const Sample& sample_2) -> distance_return_type
+{
+	PMQ pmq_1(sample_1);
+	PMQ pmq_2(sample_2);
+		
+	// find the most difference between distributions
+
+	Sample concat_data;
+	for (int i = 0; i < sample_1.size(); i++)
+	{
+		concat_data.push_back(sample_1[i]);
+	}
+	for (int i = 0; i < sample_2.size(); i++)
+	{
+		concat_data.push_back(sample_2[i]);
+	}
+	std::sort(concat_data.begin(), concat_data.end());
+	
+	D max_difference = 0;
+	for (int i = 0; i < concat_data.size(); i++)
+	{
+		if (abs(pmq_1.cdf(concat_data[i]) - pmq_2.cdf(concat_data[i])) > max_difference)
+		{
+			max_difference = abs(pmq_1.cdf(concat_data[i]) - pmq_2.cdf(concat_data[i]));
+		}
+	}
+
+
+    return max_difference;
+}
+
+
+}  // namespace metric
+#endif
diff --git a/modules/distance/k-random/KolmogorovSmirnov.hpp b/modules/distance/k-random/KolmogorovSmirnov.hpp
new file mode 100644
index 00000000..40db0302
--- /dev/null
+++ b/modules/distance/k-random/KolmogorovSmirnov.hpp
@@ -0,0 +1,50 @@
+/*
+This Source Code Form is subject to the terms of the Mozilla Public
+License, v. 2.0. If a copy of the MPL was not distributed with this
+file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+Copyright (c) 2020 Panda Team
+*/
+
+#ifndef _METRIC_DISTANCE_K_RANDOM_KOLMOGOROV_SMIRNOV_HPP
+#define _METRIC_DISTANCE_K_RANDOM_KOLMOGOROV_SMIRNOV_HPP
+
+namespace metric {
+	
+
+/**
+ * @brief
+ *
+ * To compare the two samples, we construct their empirical cumulative distribution functions (CDF). CDFs is calculated by PMQ. 
+ * The Kolmogorov-Smirnov (KS) distance is defined to be the largest absolute difference between the two empirical CDFs evaluated at any point.
+ * 
+ * @tparam Sample - sample type
+ * @tparam D - distance return type
+ */
+template <typename Sample, typename D = double>
+struct KolmogorovSmirnov {
+    using distance_return_type = D;
+
+    /**
+     * @brief Construct a new Kolmogorov-Smirnov object
+     *
+     */
+    explicit KolmogorovSmirnov()
+    {
+    }
+
+    /**
+     * @brief calculate Kolmogorov-Smirnov distance between two samples
+     *
+     * @param sample_1 first sample
+     * @param sample_2 second sample
+     * @return distance
+     */
+    distance_return_type operator()(const Sample& sample_1, const Sample& sample_2);
+};
+
+
+}  // namespace metric
+
+#include "KolmogorovSmirnov.cpp"
+#endif
diff --git a/modules/distance/k-random/RandomEMD.cpp b/modules/distance/k-random/RandomEMD.cpp
new file mode 100644
index 00000000..987ff3c0
--- /dev/null
+++ b/modules/distance/k-random/RandomEMD.cpp
@@ -0,0 +1,51 @@
+/*
+This Source Code Form is subject to the terms of the Mozilla Public
+License, v. 2.0. If a copy of the MPL was not distributed with this
+file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+Copyright (c) 2020 Panda Team
+*/
+#ifndef _METRIC_DISTANCE_K_RANDOM_EMD_CPP
+#define _METRIC_DISTANCE_K_RANDOM_EMD_CPP
+
+#include "RandomEMD.hpp"
+#include "../../utils/poor_mans_quantum.hpp"
+
+namespace metric {
+
+template <typename Sample, typename D>
+auto RandomEMD<typename Sample, typename D>::operator()(const Sample& sample_1, const Sample& sample_2) -> distance_return_type
+{
+
+	PMQ pmq_1(sample_1);
+	PMQ pmq_2(sample_2);
+		
+	// find the area between distributions
+
+	Sample concat_data;
+	for (int i = 0; i < sample_1.size(); i++)
+	{
+		concat_data.push_back(sample_1[i]);
+	}
+	for (int i = 0; i < sample_2.size(); i++)
+	{
+		concat_data.push_back(sample_2[i]);
+	}
+	std::sort(concat_data.begin(), concat_data.end());
+
+	//
+	
+	D area = 0;
+	double step = (concat_data[concat_data.size() - 1] - concat_data[0]) * precision_;
+
+	for (double value = concat_data[0]; value <= concat_data[concat_data.size() - 1]; value += step)
+	{
+		area += abs(pmq_1.cdf(value) - pmq_2.cdf(value)) * step;
+	}
+
+    return area;
+}
+
+
+}  // namespace metric
+#endif
diff --git a/modules/distance/k-random/RandomEMD.hpp b/modules/distance/k-random/RandomEMD.hpp
new file mode 100644
index 00000000..77d9f340
--- /dev/null
+++ b/modules/distance/k-random/RandomEMD.hpp
@@ -0,0 +1,64 @@
+/*
+This Source Code Form is subject to the terms of the Mozilla Public
+License, v. 2.0. If a copy of the MPL was not distributed with this
+file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+Copyright (c) 2020 Panda Team
+*/
+
+#ifndef _METRIC_DISTANCE_K_RANDOM_EMD_HPP
+#define _METRIC_DISTANCE_K_RANDOM_EMD_HPP
+
+namespace metric {
+	
+
+/**
+ * @brief
+ *
+ * Earth Mover�s Distance (EMD), also known as the first Wasserstein distance. 
+ * �hysical interpretation �� the EMDis easy to understand: imagine the two datasets to be piles of earth, and the goal is to move the first pile around to match the second. 
+ * The Earth Mover�s Distance is the minimum amount of work involved, where �amount of work� is the amount of earth you have to move multiplied by the distance you have to move it. 
+ * The EMD can also be shown to be equal to the area between the two empirical CDFs, which is calculated py PMQ. 
+ *
+ * @tparam Sample - sample type
+ * @tparam D - distance return type
+ */
+template <typename Sample, typename D = double>
+struct RandomEMD {
+    using distance_return_type = D;
+
+    /**
+     * @brief Construct a new EMD object
+     *
+     */
+    explicit RandomEMD()
+    {
+    }
+
+    /**
+     * @brief Construct a new EMD object
+     *
+     * @param precision used for integration. Should be in (0, 1). Less means more accurate.
+     */
+    explicit RandomEMD(double precision) : precision_(precision)
+    {
+    }
+
+    /**
+     * @brief calculate EMD distance between two samples
+     *
+     * @param sample_1 first sample
+     * @param sample_2 second sample
+     * @return distance
+     */
+    distance_return_type operator()(const Sample& sample_1, const Sample& sample_2);
+
+private:
+	double precision_ = 0.01;
+};
+
+
+}  // namespace metric
+
+#include "RandomEMD.cpp"
+#endif
diff --git a/modules/mapping/KOC.hpp b/modules/mapping/KOC.hpp
index 09d42e0e..92693c9c 100644
--- a/modules/mapping/KOC.hpp
+++ b/modules/mapping/KOC.hpp
@@ -51,13 +51,13 @@ namespace metric {
 			/**
 				 * @brief Construct a new KOC object
 				 * 
-				 * @param graph
-				 * @param metric 
+				 * @param graph - pre created graph (with metric::Graph interface) over which SOM is being constructed. Default is metric::Grid6 (hexagones grid).
+				 * @param metric - metric or distance being used for SOM training. Default is metric::Euclidian.
 				 * @param anomaly_sigma
-				 * @param start_learn_rate
-				 * @param finish_learn_rate
-				 * @param iterations
-				 * @param distribution
+				 * @param s_learn_rate - start learning rate for SOM training. Shows how error can influence on the sample on the first iterations of the training.
+				 * @param f_learn_rate - finish learning rate for SOM training. Shows how error can influence on the sample on the last iterations of the training.
+				 * @param iterations - maximum number of iterations for SOM training.
+				 * @param distribution - distribution used for creating initial weights (positions) of the SOM nodes. 
 				 */
 			KOC(Graph graph, Metric metric, double anomaly_sigma = 1.0, double start_learn_rate = 0.8, double finish_learn_rate = 0.0, size_t iterations = 20, Distribution distribution = Distribution(-1, 1))
 				: som_(graph, metric, start_learn_rate, finish_learn_rate, iterations, distribution),
@@ -68,16 +68,16 @@ namespace metric {
 			/**
 				 * @brief Construct a new KOC object
 				 * 
-				 * @param graph 
-				 * @param metric 
+				 * @param graph - pre created graph (with metric::Graph interface) over which SOM is being constructed. Default is metric::Grid6 (hexagones grid).
+				 * @param metric - metric or distance being used for SOM training. Default is metric::Euclidian.
 				 * @param anomaly_sigma
-				 * @param start_learn_rate
-				 * @param finish_learn_rate
-				 * @param iterations 
-				 * @param distribution 
-				 * @param neighborhood_start_size 
-				 * @param neigbour_range_decay 
-				 * @param random_seed 
+				 * @param s_learn_rate - start learning rate for SOM training. Shows how error can influence on the sample on the first iterations of the training.
+				 * @param f_learn_rate - finish learning rate for SOM training. Shows how error can influence on the sample on the last iterations of the training.
+				 * @param iterations - maximum number of iterations for SOM training.
+				 * @param distribution - distribution used for creating initial weights (positions) of the SOM nodes. 
+				 * @param neighborhood_start_size - number of neighbours of the checking node which weights (positions) will be corrected while SOM training. 
+				 * @param neigbour_range_decay - shows how distance from the checking node influence to weights (positions) of the checking node neighbours while SOM training. 
+				 * @param random_seed - the value used to initialize the random number generator, which is used for the initial distribution of weight (positions).
 				 */
 			KOC(Graph graph, Metric metric, double anomaly_sigma, double start_learn_rate, double finish_learn_rate, size_t iterations,
 				Distribution distribution, double neighborhood_start_size, double neigbour_range_decay, long long random_seed)
@@ -205,72 +205,71 @@ namespace metric {
 		typedef typename recType::value_type T;
 
 		/**
-			* @brief
-			*
-			* @param nodesNumber
-			* @param start_learn_rate
-			* @param finish_learn_rate
-			* @param iterations
-			* @param distribution_min
-			* @param distribution_max
+		 * @brief
+		 * @param nodesNumber - total number of the nodes for graph over which SOM is being constructed. From this number it should be possible to extract the int square root.
+		 * @param start_learn_rate - start learning rate for SOM training. Shows how error can influence on the sample on the first iterations of the training.
+		 * @param finish_learn_rate - finish learning rate for SOM training. Shows how error can influence on the sample on the last iterations of the training.
+		 * @param iterations - maximum number of iterations for SOM training.
+		 * @param distribution_min - left side (minimum) of the distribution used for creating initial weights (positions) of the SOM nodes. 
+		 * @param distribution_max - right side (maximum) of the distribution used for creating initial weights (positions) of the SOM nodes. 
 		 */
 		KOC_factory(size_t nodesNumber, double anomaly_sigma, double start_learn_rate = 0.8, double finish_learn_rate = 0.0, size_t iterations = 20, T distribution_min = -1, T distribution_max = 1);
 
 		/**
-			* @brief
-			*
-			* @param nodesWidth
-			* @param nodesHeight
-			* @param start_learn_rate
-			* @param finish_learn_rate
-			* @param iterations
-			* @param distribution_min
-			* @param distribution_max
+		 * @brief
+		 *
+		 * @param nodesWidth - number of the nodes in width for graph over which SOM is being constructed.
+		 * @param nodesHeight - number of the nodes in height for graph over which SOM is being constructed.
+		 * @param start_learn_rate - start learning rate for SOM training. Shows how error can influence on the sample on the first iterations of the training.
+		 * @param finish_learn_rate - finish learning rate for SOM training. Shows how error can influence on the sample on the last iterations of the training.
+		 * @param iterations - maximum number of iterations for SOM training.
+		 * @param distribution_min - left side (minimum) of the distribution used for creating initial weights (positions) of the SOM nodes. 
+		 * @param distribution_max - right side (maximum) of the distribution used for creating initial weights (positions) of the SOM nodes. 
 		 */
 		KOC_factory(size_t nodesWidth = 5, size_t nodesHeight = 4, double anomaly_sigma = 1.0, double start_learn_rate = 0.8, double finish_learn_rate = 0.0, size_t iterations = 20, T distribution_min = -1, T distribution_max = 1);
 
 		/**
-			* @brief
-			*
-			* @param nodesNumber
-			* @param start_learn_rate
-			* @param finish_learn_rate
-			* @param iterations
-			* @param distribution_min
-			* @param distribution_max
-			* @param neighborhood_start_size
-			* @param neigbour_range_decay
-			* @param random_seed
+		 * @brief
+		 *
+		 * @param nodesNumber - total number of the nodes for graph over which SOM is being constructed. From this number it should be possible to extract the int square root.
+		 * @param start_learn_rate - start learning rate for SOM training. Shows how error can influence on the sample on the first iterations of the training.
+		 * @param finish_learn_rate - finish learning rate for SOM training. Shows how error can influence on the sample on the last iterations of the training.
+		 * @param iterations - maximum number of iterations for SOM training.
+		 * @param distribution_min - left side (minimum) of the distribution used for creating initial weights (positions) of the SOM nodes. 
+		 * @param distribution_max - right side (maximum) of the distribution used for creating initial weights (positions) of the SOM nodes. 
+		 * @param neighborhood_start_size - number of neighbours of the checking node which weights (positions) will be corrected while SOM training. 
+		 * @param neigbour_range_decay - shows how distance from the checking node influence to weights (positions) of the checking node neighbours while SOM training. 
+		 * @param random_seed - the value used to initialize the random number generator, which is used for the initial distribution of weight (positions).
 		 */
 		KOC_factory(size_t nodesNumber, double anomaly_sigma,
 			double start_learn_rate, double finish_learn_rate, size_t iterations, T distribution_min, T distribution_max, 
 			double neighborhood_start_size, double neigbour_range_decay, long long random_seed);
 
 		/**
-			* @brief
-			*
-			* @param nodesWidth
-			* @param nodesHeight
-			* @param start_learn_rate
-			* @param finish_learn_rate
-			* @param iterations
-			* @param distribution_min
-			* @param distribution_max
-			* @param neighborhood_start_size
-			* @param neigbour_range_decay
-			* @param random_seed
+		 * @brief
+		 *
+		 * @param nodesWidth - number of the nodes in width for graph over which SOM is being constructed.
+		 * @param nodesHeight - number of the nodes in height for graph over which SOM is being constructed.
+		 * @param s_learn_rate - start learning rate for SOM training. Shows how error can influence on the sample on the first iterations of the training.
+		 * @param f_learn_rate - finish learning rate for SOM training. Shows how error can influence on the sample on the last iterations of the training.
+		 * @param iterations - maximum number of iterations for SOM training.
+		 * @param distribution_min - left side (minimum) of the distribution used for creating initial weights (positions) of the SOM nodes. 
+		 * @param distribution_max - right side (maximum) of the distribution used for creating initial weights (positions) of the SOM nodes. 
+		 * @param neighborhood_start_size - number of neighbours of the checking node which weights (positions) will be corrected while SOM training. 
+		 * @param neigbour_range_decay - shows how distance from the checking node influence to weights (positions) of the checking node neighbours while SOM training. 
+		 * @param random_seed - the value used to initialize the random number generator, which is used for the initial distribution of weight (positions).
 		 */
 		KOC_factory(size_t nodesWidth, size_t nodesHeight, double anomaly_sigma,
 			double start_learn_rate, double finish_learn_rate, size_t iterations, T distribution_min, T distribution_max, 
 			double neighborhood_start_size, double neigbour_range_decay, long long random_seed);
 
 		/**
-			* @brief
-			*
-			* @param samples
-			* @param num_clusters
-			*
-			* @return KOC_details::KOC<recType, Graph, Metric, Distribution>
+		 * @brief 
+		 * 
+		 * @param samples
+		 * @param num_clusters - number of expected clusters
+		 *
+		 * @return KOC_details::KOC<recType, Graph, Metric, Distribution>
 		 */
 		KOC_details::KOC<recType, Graph, Metric, Distribution> operator()(const std::vector<recType>& samples, int num_clusters, int min_cluster_size = 1);
 
diff --git a/modules/utils/poor_mans_quantum.hpp b/modules/utils/poor_mans_quantum.hpp
index 9e3727f8..37f89485 100644
--- a/modules/utils/poor_mans_quantum.hpp
+++ b/modules/utils/poor_mans_quantum.hpp
@@ -47,22 +47,24 @@ class PMQ {
      * @param d 
      */
     PMQ(T par1 = 1.0, T par2 = 0.0, size_t n = 1000, Distribution d = Distribution());
-
+	
+	// Updated by Stepan Mamontov 26.02.2020
     /**
      * @brief Construct a new PMQ object
      * 
      * @param data 
      * @param d 
      */
-    PMQ(std::vector<double> data, Distribution d = Distribution());
-
+    PMQ(std::vector<T> data, Distribution d = Distribution());
+	
+	// Updated by Stepan Mamontov 26.02.2020
     /**
      * @brief Construct a new PMQ object
      * 
      * @param data 
      * @param d 
      */
-    PMQ(std::vector<float> data, Distribution d = Distribution());
+    //PMQ(std::vector<float> data, Distribution d = Distribution());
 
     /**
      * @brief 
diff --git a/modules/utils/poor_mans_quantum/distributions/Discrete.cpp b/modules/utils/poor_mans_quantum/distributions/Discrete.cpp
index 549ab3f3..f4ad64cf 100644
--- a/modules/utils/poor_mans_quantum/distributions/Discrete.cpp
+++ b/modules/utils/poor_mans_quantum/distributions/Discrete.cpp
@@ -49,22 +49,25 @@ T Discrete<T>::quantil(T p)
 template <typename T>
 T Discrete<T>::mean()
 {
-    T sum = 0;
+	// Updated mean calculation by Stepan Mamontov 04.02.2020
+	// Updated according to formula: mean = sum (x * P(x))  ->  https://en.wikipedia.org/wiki/Mean#Mean_of_a_probability_distribution
+    T sum = _data[0] * _prob[0];
     for (size_t i = 1; i < _data.size(); ++i) {
-        sum += (_data[i] - _data[i - 1]) * (_prob[i] - _prob[i - 1]);
+        sum += _data[i] * (_prob[i] - _prob[i - 1]);
     }
-    T value = sum * (_data[_data.size() - 1] - _data[0]);  // / T(_data.size());
 
-    return value;
+    return sum;
 }
 
 template <typename T>
 T Discrete<T>::variance()
 {
+	// Updated variance calculation by Stepan Mamontov 04.02.2020
+	// Updated accroding to formula: variance = sum ((x - mean) * (x - mean) * P(x))  ->  https://en.wikipedia.org/wiki/Variance#Discrete_random_variable
     T mean_value = mean();
-    T sum = 0;
-    for (size_t i = 0; i < _data.size(); ++i) {
-        sum += (_data[i] - mean_value) * (_data[i] - mean_value);
+    T sum = (_data[0] - mean_value) * (_data[0] - mean_value) * _prob[0];
+    for (size_t i = 1; i < _data.size(); ++i) {
+        sum += (_data[i] - mean_value) * (_data[i] - mean_value) * (_prob[i] - _prob[i-1]);
     }
     return sum;
 }
@@ -73,7 +76,24 @@ T Discrete<T>::variance()
 template <typename T>
 T Discrete<T>::cdf(const T x)
 {
-    return (T(0.5) * (T(1) + (x - _p1) / (_p2 * 1.41421356237309504880)));
+	// Updated cdf calculation by Stepan Mamontov 04.02.2020
+	// Updated by return interpolated func of cumulative propability from distributed value
+    std::vector<T> values;
+    values.push_back(x);
+	T result = akimaInterp1(_data, _prob, values)[0];
+	
+	// cut probs over the _data values
+	if (x < _data[0])
+	{
+		result = 0;
+	}
+	
+	if (x > _data[_data.size() - 1])
+	{
+		result = 1;
+	}
+
+	return result;
 }
 
 }  // end namespace metric
diff --git a/modules/utils/poor_mans_quantum/poor_mans_quantum.cpp b/modules/utils/poor_mans_quantum/poor_mans_quantum.cpp
index 1d0b3283..c13ca4cc 100644
--- a/modules/utils/poor_mans_quantum/poor_mans_quantum.cpp
+++ b/modules/utils/poor_mans_quantum/poor_mans_quantum.cpp
@@ -123,37 +123,58 @@ PMQ<Distribution, T>::PMQ(T par1, T par2, size_t n, Distribution d)
     }
 }
 
+// Updated by Stepan Mamontov 26.02.2020 (updated to use generic types)
 /*** constructor for discrete samples ***/
-template <typename Distribution, typename T>
-PMQ<Distribution, T>::PMQ(std::vector<float> data, Distribution d)
-    : _dist(d)
-    , _generator(std::random_device {}())
-{
-    _dist._data = data;
-#if USE_VECTOR_SORT
-    vector_sort::sort(_dist._data);
-#else
-    std::sort(_dist._data.begin(), _dist._data.end());
-#endif
-    _dist._prob = linspace(T(0.5) / T(data.size()), T(1) - T(0.5) / T(data.size()), data.size());
-}
-
+//template <typename Distribution, typename T>
+//PMQ<Distribution, T>::PMQ(std::vector<T> data, Distribution d)
+//    : _dist(d)
+//    , _generator(std::random_device {}())
+//{
+//    _dist._data.resize(data.size());
+//    for (size_t i = 0; i < data.size(); ++i) {
+//        _dist._data[i] = T(data[i]);
+//    }
+//#if USE_VECTOR_SORT
+//    vector_sort::sort(_dist._data);
+//#else
+//    std::sort(_dist._data.begin(), _dist._data.end());
+//#endif
+//	auto prob = linspace(T(0.5) / T(data.size()), T(1) - T(0.5) / T(data.size()), data.size());
+//
+//    _dist._prob.resize(prob.size());
+//    for (size_t i = 0; i < prob.size(); ++i) {
+//        _dist._prob[i] = T(prob[i]);
+//    }
+//}
+
+// Updated by Stepan Mamontov 26.02.2020 (calculate histogram of the data and then calculate cumulative distribution)
 /*** constructor for discrete samples ***/
 template <typename Distribution, typename T>
-PMQ<Distribution, T>::PMQ(std::vector<double> data, Distribution d)
+PMQ<Distribution, T>::PMQ(std::vector<T> data, Distribution d)
     : _dist(d)
     , _generator(std::random_device {}())
 {
-    _dist._data.resize(data.size());
-    for (size_t i = 0; i < data.size(); ++i) {
-        _dist._data[i] = float(data[i]);
-    }
 #if USE_VECTOR_SORT
-    vector_sort::sort(_dist._data);
+	vector_sort::sort(data);
 #else
-    std::sort(_dist._data.begin(), _dist._data.end());
+	std::sort(data.begin(), data.end());
 #endif
-    _dist._prob = linspace(T(0.5) / T(data.size()), T(1) - T(0.5) / T(data.size()), data.size());
+
+	std::map<T, int> hist{};
+	for (int i = 0; i < data.size(); i++)
+	{
+        ++hist[data[i]];
+    }
+	
+	_dist._data.clear();
+	_dist._prob.clear();
+	T cumulative_sum = 0;
+	for(std::map<T, int>::iterator it = hist.begin(); it != hist.end(); ++it) 
+	{
+		_dist._data.push_back(it->first);
+		cumulative_sum += (T) it->second / data.size();
+		_dist._prob.push_back(cumulative_sum);
+	}
 }
 
 template <typename Distribution, typename T>
diff --git a/tests/distance_tests/CMakeLists.txt b/tests/distance_tests/CMakeLists.txt
index a9904423..33c477ff 100644
--- a/tests/distance_tests/CMakeLists.txt
+++ b/tests/distance_tests/CMakeLists.txt
@@ -1,31 +1,61 @@
+cmake_minimum_required(VERSION 3.10)
 set(BOOST_TEST_LIB ${Boost_LIBRARIES})
 
+enable_testing()
+
+set(Boost_DEBUG on)  
+find_package(Boost)
+
+project(distance_tests)
+
+include_directories(
+	${PROJECT_SOURCE_DIR} 
+	${PROJECT_SOURCE_DIR}/../..
+    ${Boost_INCLUDE_DIRS}
+	)	
+
+# for Visual Studio
+if(MSVC) 
+	set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} /std:c++17")  
+endif(MSVC)
+
+set(BOOST_LINK_DIRS "${BOOST_ROOT}/stage/lib")
+
+# file(COPY "assets" DESTINATION ".") 
+
 if(UNIX)
-  set(THREADS_PREFER_PTHREAD_FLAG ON)
-  find_package(Threads REQUIRED)
+    set(THREADS_PREFER_PTHREAD_FLAG ON)
+    find_package(Threads REQUIRED)
+
 endif(UNIX)
 
 AUX_SOURCE_DIRECTORY(. TEST_SRCS)
 #Run through each source
 foreach(testSrc ${TEST_SRCS})
 
-  #Extract the filename without an extension (NAME_WE)
-  get_filename_component(testName ${testSrc} NAME_WE)
-
-  #Add compile target
-  add_executable(${testName} ${testSrc})
-
-  #link to Boost libraries AND your targets and dependencies
-  if(WIN32)
-    target_compile_options(${testName} PRIVATE /D BOOST_ALL_NO_LIB=1)
-  endif()
-  target_link_libraries(${testName} ${BOOST_TEST_LIB})
-
-  if(UNIX)
-    target_link_libraries(${testName} Threads::Threads)
-# uncomment next two lines if want to enable address sanitizer in debug mode
-#    set (CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -fno-omit-frame-pointer -fsanitize=address")
-#    set (CMAKE_LINKER_FLAGS_DEBUG "${CMAKE_LINKER_FLAGS_DEBUG} -fno-omit-frame-pointer -fsanitize=address")
-  endif(UNIX)
-  add_test(NAME ${testName}  COMMAND ${testName} )
+        #Extract the filename without an extension (NAME_WE)
+        get_filename_component(testName ${testSrc} NAME_WE)
+
+        #Add compile target
+        add_executable(${testName} ${testSrc})
+ 
+		#link to Boost libraries AND your targets and dependencies 
+		target_link_directories(${testName} PRIVATE ${BOOST_LINK_DIRS})
+		target_link_libraries(${testName} ${Boost_LIBRARIES})
+        target_link_libraries(${testName} ${BOOST_TEST_LIB})
+		set_target_properties(${testName} PROPERTIES CXX_STANDARD 17)
+	if(UNIX)
+		set(CMAKE_CXX_FLAGS "-O3")
+
+		target_link_libraries(${exampleName} Threads::Threads)
+
+		target_compile_options(${exampleName} PRIVATE -pthread)
+		set_target_properties(${exampleName} PROPERTIES LINK_FLAGS -pthread)
+	endif(UNIX)
+        #I like to move testing binaries into a testBin directory
+#        set_target_properties(${testName} PROPERTIES RUNTIME_OUTPUT_DIRECTORY  ${CMAKE_CURRENT_SOURCE_DIR})
+        #Finally add it to test execution - 
+        #Notice the WORKING_DIRECTORY and COMMAND
+        add_test(NAME ${testName} COMMAND ${testName})
+
 endforeach(testSrc)
diff --git a/tests/distance_tests/README.md b/tests/distance_tests/README.md
new file mode 100644
index 00000000..1589447d
--- /dev/null
+++ b/tests/distance_tests/README.md
@@ -0,0 +1,30 @@
+
+---
+
+## Run
+*You need STL and C++17 support to compile.*
+
+#### Using CMake
+
+_Windows_
+
+```bash
+mkdir build
+cd build
+cmake .. -A x64 -DBoost_NAMESPACE="libboost" -DBoost_COMPILER="-vc141" -DBOOST_ROOT=C:/boost_1_70_0/
+```
+
+Where `DBoost_NAMESPACE` and `DBoost_COMPILER` params point to suffixes and postfixes of the search libs in the `stage/lib` folder of the boost root directory. 
+And `DBOOST_ROOT` point to the boost root folder. 
+
+Then open solution in the Microsoft Visual Studio.
+
+_Linux_
+
+Just run cmake
+```bash
+mkdir build
+cd build
+cmake ..
+make
+```
\ No newline at end of file
diff --git a/tests/distance_tests/cramervon_mises_tests.cpp b/tests/distance_tests/cramervon_mises_tests.cpp
new file mode 100644
index 00000000..3e1832af
--- /dev/null
+++ b/tests/distance_tests/cramervon_mises_tests.cpp
@@ -0,0 +1,92 @@
+/*
+This Source Code Form is subject to the terms of the Mozilla Public
+License, v. 2.0. If a copy of the MPL was not distributed with this
+file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+Copyright (c) 2020 Panda Team
+*/
+
+#include <vector>
+#include "modules/distance.hpp"
+
+#define BOOST_TEST_MODULE Main
+
+#include <boost/test/unit_test.hpp>
+
+BOOST_AUTO_TEST_CASE(basic_use) 
+{
+	/*** here are some data records ***/
+	std::vector<double> samples_1 = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 };
+	std::vector<double> samples_2 = { 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11 };
+	   
+    metric::CramervonNises<std::vector<double>, double> distance;
+
+    auto result = distance(samples_1, samples_2);
+	
+    std::cout << "basic use result: " << result << std::endl;
+
+    double t = 5.0; // 5.0%
+    BOOST_CHECK_CLOSE(result, 0.311, t);
+}
+
+BOOST_AUTO_TEST_CASE(different_dimensions) 
+{
+	/*** here are some data records ***/
+	std::vector<double> samples_1 = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 };
+	std::vector<double> samples_2 = { 1, 2, 3, 4, 5 };
+	   
+    metric::CramervonNises<std::vector<double>, double> distance;
+
+    auto result = distance(samples_1, samples_2);
+	
+    std::cout << "different dimensions result: " << result << std::endl;
+	
+    double t = 5.0; // 5.0%
+    BOOST_CHECK_CLOSE(result, 0.81, t);
+}
+
+BOOST_AUTO_TEST_CASE(equal_samples) 
+{
+	/*** here are some data records ***/
+	std::vector<double> samples_1 = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 };
+	std::vector<double> samples_2 = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 };
+	   
+    metric::CramervonNises<std::vector<double>, double> distance;
+
+    auto result = distance(samples_1, samples_2);
+	
+    std::cout << "equal samples result: " << result << std::endl;
+	
+    double t = 0.01; // 0.01%
+    BOOST_CHECK_CLOSE(result, 0.0, t);
+}
+
+BOOST_AUTO_TEST_CASE(single_value_distribution) 
+{
+	/*** here are some data records ***/
+	std::vector<double> samples_1 = { 0, 0, 0, 0, 0, 0, 0 };
+	std::vector<double> samples_2 = { 1, 1, 1, 1, 1, 1, 1 };
+	   
+    metric::CramervonNises<std::vector<double>, double> distance;
+
+    auto result = distance(samples_1, samples_2);
+	
+    std::cout << "single value distribution result: " << result << std::endl;
+
+    BOOST_CHECK( std::isnan( result ) );
+}
+
+BOOST_AUTO_TEST_CASE(non_intersect_distribution) 
+{
+	/*** here are some data records ***/
+	std::vector<double> samples_1 = { 0, 1, 1, 2, 2, 2, 2, 1, 1, 0 };
+	std::vector<double> samples_2 = { 10, 11, 11, 12, 12, 12, 12, 11, 11, 10 };
+	   
+    metric::CramervonNises<std::vector<double>, double> distance;
+
+    auto result = distance(samples_1, samples_2);
+	
+    std::cout << "non intersect distributionn result: " << result << std::endl;
+
+    BOOST_CHECK( std::isnan( result ) );
+}
diff --git a/tests/distance_tests/kolmogorov_smirnov_tests.cpp b/tests/distance_tests/kolmogorov_smirnov_tests.cpp
new file mode 100644
index 00000000..72ca66ae
--- /dev/null
+++ b/tests/distance_tests/kolmogorov_smirnov_tests.cpp
@@ -0,0 +1,95 @@
+/*
+This Source Code Form is subject to the terms of the Mozilla Public
+License, v. 2.0. If a copy of the MPL was not distributed with this
+file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+Copyright (c) 2020 Panda Team
+*/
+
+#include <vector>
+#include "modules/distance.hpp"
+
+#define BOOST_TEST_MODULE Main
+
+#include <boost/test/unit_test.hpp>
+
+BOOST_AUTO_TEST_CASE(basic_use) 
+{
+	/*** here are some data records ***/
+	std::vector<double> samples_1 = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 };
+	std::vector<double> samples_2 = { 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11 };
+	   
+    metric::KolmogorovSmirnov<std::vector<double>, double> distance;
+
+    auto result = distance(samples_1, samples_2);
+	
+    std::cout << "basic use result: " << result << std::endl;
+
+    double t = 10.0; // 10.0%
+    BOOST_CHECK_CLOSE(result, 0.1, t);
+}
+
+BOOST_AUTO_TEST_CASE(different_dimensions) 
+{
+	/*** here are some data records ***/
+	std::vector<double> samples_1 = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 };
+	std::vector<double> samples_2 = { 1, 2, 3, 4, 5 };
+	   
+    metric::KolmogorovSmirnov<std::vector<double>, double> distance;
+
+    auto result = distance(samples_1, samples_2);
+	
+    std::cout << "different dimensions result: " << result << std::endl;
+
+    double t = 10.0; // 10.0%
+    BOOST_CHECK_CLOSE(result, 0.45, t);
+}
+
+BOOST_AUTO_TEST_CASE(equal_samples) 
+{
+	/*** here are some data records ***/
+	std::vector<double> samples_1 = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 };
+	std::vector<double> samples_2 = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 };
+	   
+    metric::KolmogorovSmirnov<std::vector<double>, double> distance;
+
+    auto result = distance(samples_1, samples_2);
+	
+    std::cout << "equal samples result: " << result << std::endl;
+
+    double t = 0.01; // 0.01%
+    BOOST_CHECK_CLOSE(result, 0.0, t);
+}
+
+BOOST_AUTO_TEST_CASE(single_value_distribution) 
+{
+	/*** here are some data records ***/
+	std::vector<double> samples_1 = { 0, 0, 0, 0, 0, 0, 0 };
+	std::vector<double> samples_2 = { 1, 1, 1, 1, 1, 1, 1 };
+	   
+    metric::KolmogorovSmirnov<std::vector<double>, double> distance;
+
+    auto result = distance(samples_1, samples_2);
+	
+    std::cout << "single value distribution result: " << result << std::endl;
+
+    double t = 0.01; // 0.01%
+    BOOST_CHECK_CLOSE(result, 0.0, t);
+}
+
+
+BOOST_AUTO_TEST_CASE(non_intersect_distribution) 
+{
+	/*** here are some data records ***/
+	std::vector<double> samples_1 = { 0, 1, 1, 2, 2, 2, 2, 1, 1, 0 };
+	std::vector<double> samples_2 = { 10, 11, 11, 12, 12, 12, 12, 11, 11, 10 };
+	   
+    metric::KolmogorovSmirnov<std::vector<double>, double> distance;
+
+    auto result = distance(samples_1, samples_2);
+	
+    std::cout << "non intersect distributionn result: " << result << std::endl;
+	
+    double t = 0.01; // 0.01%
+    BOOST_CHECK_CLOSE(result, 0.0, t);
+}
\ No newline at end of file
diff --git a/tests/distance_tests/random_emd_tests.cpp b/tests/distance_tests/random_emd_tests.cpp
new file mode 100644
index 00000000..518bfbaa
--- /dev/null
+++ b/tests/distance_tests/random_emd_tests.cpp
@@ -0,0 +1,91 @@
+/*
+This Source Code Form is subject to the terms of the Mozilla Public
+License, v. 2.0. If a copy of the MPL was not distributed with this
+file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+Copyright (c) 2020 Panda Team
+*/
+
+#include <vector>
+#include "modules/distance.hpp"
+
+#define BOOST_TEST_MODULE Main
+
+#include <boost/test/unit_test.hpp>
+
+BOOST_AUTO_TEST_CASE(basic_use) 
+{
+	/*** here are some data records ***/
+	std::vector<double> samples_1 = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 };
+	std::vector<double> samples_2 = { 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11 };
+	   
+    metric::RandomEMD<std::vector<double>, double> distance;
+
+    auto result = distance(samples_1, samples_2);
+	
+    std::cout << "basic use result: " << result << std::endl;
+
+    double t = 1.0; // 1.0%
+    BOOST_CHECK_CLOSE(result, 1.0, t);
+}
+BOOST_AUTO_TEST_CASE(different_dimensions) 
+{
+	/*** here are some data records ***/
+	std::vector<double> samples_1 = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 };
+	std::vector<double> samples_2 = { 1, 2, 3, 4, 5 };
+	   
+    metric::RandomEMD<std::vector<double>, double> distance;
+
+    auto result = distance(samples_1, samples_2);
+	
+    std::cout << "different dimensions result: " << result << std::endl;
+	
+    double t = 1.0; // 1.0%
+    BOOST_CHECK_CLOSE(result, 2.18, t);
+}
+
+BOOST_AUTO_TEST_CASE(equal_samples) 
+{
+	/*** here are some data records ***/
+	std::vector<double> samples_1 = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 };
+	std::vector<double> samples_2 = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 };
+	   
+    metric::RandomEMD<std::vector<double>, double> distance;
+
+    auto result = distance(samples_1, samples_2);
+	
+    std::cout << "equal samples result: " << result << std::endl;
+	
+    double t = 0.01; // 0.01%
+    BOOST_CHECK_CLOSE(result, 0.0, t);
+}
+
+BOOST_AUTO_TEST_CASE(single_value_distribution) 
+{
+	/*** here are some data records ***/
+	std::vector<double> samples_1 = { 0, 0, 0, 0, 0, 0, 0 };
+	std::vector<double> samples_2 = { 1, 1, 1, 1, 1, 1, 1 };
+	   
+    metric::RandomEMD<std::vector<double>, double> distance;
+
+    auto result = distance(samples_1, samples_2);
+	
+    std::cout << "single value distribution result: " << result << std::endl;
+
+    BOOST_CHECK( std::isnan( result ) );
+}
+
+BOOST_AUTO_TEST_CASE(non_intersect_distribution) 
+{
+	/*** here are some data records ***/
+	std::vector<double> samples_1 = { 0, 1, 1, 2, 2, 2, 2, 1, 1, 0 };
+	std::vector<double> samples_2 = { 10, 11, 11, 12, 12, 12, 12, 11, 11, 10 };
+	   
+    metric::RandomEMD<std::vector<double>, double> distance;
+
+    auto result = distance(samples_1, samples_2);
+	
+    std::cout << "non intersect distributionn result: " << result << std::endl;
+
+    BOOST_CHECK( std::isnan( result ) );
+}

From f732e26278ccca45197e93efedff2c206dca2049 Mon Sep 17 00:00:00 2001
From: Stepka <micromillioner@yahoo.com>
Date: Sat, 21 Mar 2020 00:26:19 +0200
Subject: [PATCH 4/4] Knn search graph and search (#81)

* [KNN graph] - added KNN graph and its construction using random pair division

* [KNN] - gnn search blanks

* [KNN] added knn graph search

* [KNN] added some comments

* [KNN] added tests for knn graph

* [KNN Graph] added copy constructor

* [KNN Graph] added constructor from Tree object

* [KNN Graph] update knn-graph constructor from tree (optimized)

* [KNN Graph] Updated tests

* [KNN Graph] little cleanup

* [KNN] lilfix

* [KNN] resolved merge conflicts
---
 examples/space_examples/CMakeLists.txt        |   6 +-
 examples/space_examples/knn_graph_example.cpp | 102 +++++
 examples/space_examples/simple_example.cpp    |   3 +
 modules/space/tree.cpp                        |  22 +
 modules/space/tree.hpp                        |   9 +
 modules/utils/graph.hpp                       |  99 +++++
 modules/utils/graph/graph.cpp                 | 386 ++++++++++++++++++
 tests/mapping_tests/CMakeLists.txt            |   9 +-
 tests/mapping_tests/mapping_tests.cpp         |   2 +-
 tests/space_tests/CMakeLists.txt              |  39 +-
 tests/space_tests/README.md                   |  33 ++
 tests/space_tests/knn_graph_test.cpp          | 313 ++++++++++++++
 12 files changed, 1008 insertions(+), 15 deletions(-)
 create mode 100644 examples/space_examples/knn_graph_example.cpp
 create mode 100644 tests/space_tests/README.md
 create mode 100644 tests/space_tests/knn_graph_test.cpp

diff --git a/examples/space_examples/CMakeLists.txt b/examples/space_examples/CMakeLists.txt
index a5c84998..2125f704 100644
--- a/examples/space_examples/CMakeLists.txt
+++ b/examples/space_examples/CMakeLists.txt
@@ -8,12 +8,12 @@ project(space_examples)
 file(GLOB EXAMPLE_SRCS RELATIVE ${CMAKE_CURRENT_SOURCE_DIR} *.cpp)
 
 
-include_directories(${Boost_INCLUDE_DIRS})
+include_directories(${Boost_INCLUDE_DIRS}) 
 
 # for Visual Studio
 if(MSVC)
 	set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} /std:c++17")
-endif(MSVC)
+endif(MSVC) 
  
 #Run through each source
 foreach(exampleSrc ${EXAMPLE_SRCS})
@@ -22,7 +22,7 @@ foreach(exampleSrc ${EXAMPLE_SRCS})
   get_filename_component(exampleName ${exampleSrc} NAME_WE) 
 
   #Add compile target
-  add_executable(${exampleName} ${exampleSrc})
+  add_executable(${exampleName} ${exampleSrc}) 
   
  
   #link to Boost libraries AND your targets and dependencies
diff --git a/examples/space_examples/knn_graph_example.cpp b/examples/space_examples/knn_graph_example.cpp
new file mode 100644
index 00000000..40527d77
--- /dev/null
+++ b/examples/space_examples/knn_graph_example.cpp
@@ -0,0 +1,102 @@
+/*
+This Source Code Form is subject to the terms of the Mozilla Public
+License, v. 2.0. If a copy of the MPL was not distributed with this
+file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+Copyright (c) 2019 Panda Team
+*/
+
+#include <iostream>
+#include "../../modules/utils/graph.hpp"
+#include "../../modules/distance.hpp"
+
+
+int main()
+{
+	std::cout << "Graph space example have started" << std::endl;
+	std::cout << "" << std::endl;
+	
+    size_t neighbors_num = 3;
+		
+    std::vector<std::vector<double>> table = {
+        { 0, 1 },
+        { 1, 1 },
+        { 2, 2 },
+        { 3, 3 },
+        { 4, 3 },
+        { 5, 3 },
+        { 4, 6 },
+        { 5, 1 },
+        { 4, 1 },
+        { 3, 2 },
+        { 0, 3 },
+        { 1, 3 },
+        { 2, 3 },
+        { 6, 6 },
+        { 0, 2 },
+        { 0, 9 },
+        { 0, 4 },
+        { 0, 5 },
+        { 0, 6 },
+        { 0, 7 },
+        { 0, 8 },
+    };
+	
+    auto g = metric::KNNGraph<std::vector<double>, metric::Euclidian<double>>(table, neighbors_num, 2.5 * neighbors_num);
+	
+    std::cout << "graph:" << std::endl;
+    std::cout << g.get_matrix() << std::endl;
+	
+	//
+	
+	std::vector<double> query = { 7, 5 };
+
+    std::cout << std::endl;
+    std::cout << "gnn search:" << std::endl;
+	auto found = g.gnnn_search(query, 3);
+	
+	for (size_t i = 0; i < found.size(); i++)
+	{
+		std::cout << found[i] << " -> ";
+		print_vector(g.get_node_data(found[i]));
+		std::cout << std::endl;
+	}
+    std::cout << std::endl;
+
+	//
+	
+    auto other_g = metric::KNNGraph(g);
+	
+    std::cout << "other graph:" << std::endl;
+    std::cout << "gnn search:" << std::endl;
+	found = other_g.gnnn_search(query, 3);
+		
+	for (size_t i = 0; i < found.size(); i++)
+	{
+		std::cout << found[i] << " -> ";
+		print_vector(other_g.get_node_data(found[i]));
+		std::cout << std::endl;
+	}
+    std::cout << std::endl;
+
+	//
+
+ //   metric::Tree<std::vector<double>, metric::Euclidian<double>> tree(table);
+	//
+ //   auto graph_from_tree = metric::KNNGraph(tree, neighbors_num, 2.5 * neighbors_num);
+	//
+ //   std::cout << "tree graph:" << std::endl;
+ //   std::cout << "gnn search:" << std::endl;
+	//found = graph_from_tree.gnnn_search(query, 3);
+	//
+	//for (size_t i = 0; i < found.size(); i++)
+	//{
+	//	std::cout << found[i] << " -> ";
+	//	print_vector(graph_from_tree.get_node_data(found[i]));
+	//	std::cout << std::endl;
+	//}
+ //   std::cout << std::endl;
+
+	
+    return 0;
+}
diff --git a/examples/space_examples/simple_example.cpp b/examples/space_examples/simple_example.cpp
index 2335e386..6506cfd2 100644
--- a/examples/space_examples/simple_example.cpp
+++ b/examples/space_examples/simple_example.cpp
@@ -42,6 +42,9 @@ int main()
 
     cTree.print();
 
+    std::cout << std::endl;
+    std::cout << std::endl;
+
 	//
 
     /*** batch insert ***/
diff --git a/modules/space/tree.cpp b/modules/space/tree.cpp
index 04b65af6..6bd6617b 100644
--- a/modules/space/tree.cpp
+++ b/modules/space/tree.cpp
@@ -895,6 +895,28 @@ void Tree<recType, Metric>::print_(NodeType* node_p, std::ostream& ostr) const
     }
 }
 
+template <class recType, class Metric>
+auto Tree<recType, Metric>::get_all_nodes() -> std::vector<Node_ptr>
+{	
+	std::vector<Node_ptr> all_nodes;
+	
+	all_nodes.push_back(root);
+	get_all_nodes_(root, all_nodes);
+
+	return all_nodes;
+}
+
+template <class recType, class Metric>
+void Tree<recType, Metric>::get_all_nodes_(Node_ptr node_p, std::vector<Node_ptr>& output)
+{	
+	for (std::size_t i = 0; i < node_p->children.size(); ++i) 
+	{
+		output.push_back(node_p->children[i]);
+		get_all_nodes_(node_p->children[i], output);
+	}
+}
+
+
 /*** traverse the tree from root and do something with every node ***/
 template <class recType, class Metric>
 void Tree<recType, Metric>::traverse(const std::function<void(Node_ptr)>& f)
diff --git a/modules/space/tree.hpp b/modules/space/tree.hpp
index 7026588e..834144d4 100644
--- a/modules/space/tree.hpp
+++ b/modules/space/tree.hpp
@@ -349,6 +349,13 @@ class Tree {
      */
     bool same_tree(const Node_ptr lhs, const Node_ptr rhs) const;
 
+    /**
+     * @brief recursively iterate through the tree and return all nodes of the tree
+     *
+     * @return return all nodes of the tree
+     */	
+	auto Tree<recType, Metric>::get_all_nodes() -> std::vector<Node_ptr>;
+
     /**
      * @brief compare tree with another
      *
@@ -491,6 +498,8 @@ class Tree {
     std::pair<Distance, std::size_t> distance_to_level(Node_ptr &p, int level) const;
     Distance distance_by_node(Node_ptr p1, Node_ptr p2) const;
     std::pair<Distance, std::size_t> graph_distance(Node_ptr p1, Node_ptr p2) const; 
+
+	void get_all_nodes_(Node_ptr node_p, std::vector<Node_ptr>& output);
 };
 }  // namespace metric
 #include "tree.cpp"  // include the implementation
diff --git a/modules/utils/graph.hpp b/modules/utils/graph.hpp
index c656d557..9446598a 100644
--- a/modules/utils/graph.hpp
+++ b/modules/utils/graph.hpp
@@ -10,6 +10,7 @@ Copyright (c) 2019 Michael Welsch
 #define _METRIC_UTILS_GRAPH_HPP
 
 #include "../../3rdparty/blaze/Blaze.h"
+#include "../space/tree.hpp"
 
 #include <stack>
 #include <type_traits>
@@ -62,6 +63,7 @@ class Graph {
      * @param edgesPairs 
      */
     Graph(const std::vector<std::pair<size_t, size_t>> & edgesPairs);
+
     /**
      * @brief Construct a new Graph object
      * 
@@ -325,6 +327,103 @@ class RandomUniform : public Graph<WType, isDense, false> {
 
 
 
+/**
+ * @class KNNGraph
+ * @brief Fast hierarchial method algorithm that constructs an approximate kNN graph. The method is simple and it works with any type of data for which a distance function can be provided. 
+ * Algorith has two parts. In the first part we create a crude approximation of the graph by subdividing the dataset until each subset reaches given max_bruteforce_size.
+ * In the second part this approximation is iteratively fine-tuned by combining the first algorithm with NN-descent method. 
+ * 
+ */
+template <typename Sample, typename Distance, typename WeightType = bool, bool isDense = false, bool isSymmetric = true>
+class KNNGraph : public Graph<WeightType, isDense, isSymmetric> {
+public:
+
+    /**
+     * @brief Construct a new KNN Graph object
+     * 
+     */
+    KNNGraph(std::vector<Sample> X, size_t neighbors_num, size_t max_bruteforce_size, int max_iterations = 100, double update_range = 0.02);
+
+    /**
+     * @brief Construct a new KNN Graph object
+     * 
+     */
+    KNNGraph(const KNNGraph& graph);
+
+    ///**
+    // * @brief Construct a new KNN Graph object
+    // * 
+    // */
+    //KNNGraph(Tree<Sample, Distance>& tree, size_t neighbors_num, size_t max_bruteforce_size, int max_iterations = 100, double update_range = 0.02);
+
+    ///**
+    // * @brief Construct a new KNN Graph object
+    // * 
+    // */
+    //KNNGraph(std::vector<std::vector<typename Distance::value_type>> distance_matrix);
+
+	
+	size_t neighbors_num() { return _neighbors_num; };
+
+	
+	Sample get_node_data(size_t i) { return _nodes[i]; };
+
+	
+	std::vector<int> gnnn_search(Sample query, int max_closest_num, int iterations = 10, int num_greedy_moves = -1, int num_expansions = -1);
+	
+protected:
+	size_t _neighbors_num = 1;
+	size_t _max_bruteforce_size = 10;
+	
+	int _max_iterations = 100;
+	double _update_range = 0.02;
+
+	bool _not_more_neighbors = false;
+
+	std::vector<Sample> _nodes;
+	std::vector<std::vector<typename Distance::value_type>> _distance_matrix;
+
+private:
+    /**
+     * @brief 
+     * 
+     */
+    void construct(std::vector<Sample> X);
+
+    /**
+     * @brief 
+     * 
+     */
+    void calculate_distance_matrix(std::vector<Sample> X);
+
+    /**
+     * @brief 
+     * 
+     */
+    void make_edge_pairs(std::vector<Sample> X);
+	
+    /**
+     * @brief 
+     * 
+     */
+	std::vector<std::pair<size_t, size_t>> brute_force(std::vector<Sample> samples, std::vector<int> ids);
+	
+    /**
+     * @brief 
+     * 
+     */
+	std::vector<std::pair<size_t, size_t>> random_pair_division(std::vector<Sample> samples, std::vector<int> ids, int max_size);
+
+    /**
+     * @brief 
+     * 
+     */
+	template <typename T1>
+	std::vector<size_t> sort_indexes(const std::vector<T1> &v);
+};
+
+
+
 
 
 /**
diff --git a/modules/utils/graph/graph.cpp b/modules/utils/graph/graph.cpp
index 4dc0c25a..819fcc9f 100644
--- a/modules/utils/graph/graph.cpp
+++ b/modules/utils/graph/graph.cpp
@@ -675,6 +675,392 @@ void RandomUniform<WType, isDense>::fill(MType& matrix, WType lower_bound, WType
         }
 }
 
+// KNN-Graph
+
+template <typename Sample, typename Distance, typename WeightType, bool isDense, bool isSymmetric>
+KNNGraph<Sample, Distance, WeightType, isDense, isSymmetric>::KNNGraph(std::vector<Sample> samples, size_t neighbors_num, size_t max_bruteforce_size, int max_iterations, double update_range)
+    : Graph<WeightType, isDense, isSymmetric>(samples.size()), _nodes(samples), _neighbors_num(neighbors_num), _max_bruteforce_size(max_bruteforce_size), _max_iterations(max_iterations), _update_range(update_range)
+{
+    construct(samples);
+}
+
+template <typename Sample, typename Distance, typename WeightType, bool isDense, bool isSymmetric>
+KNNGraph<Sample, Distance, WeightType, isDense, isSymmetric>::KNNGraph(const KNNGraph& graph)
+    : Graph<WeightType, isDense, isSymmetric>(), _neighbors_num(graph._neighbors_num), _max_bruteforce_size(graph._max_bruteforce_size), _max_iterations(graph._max_iterations), _update_range(graph._update_range)
+{
+	_nodes = graph._nodes;
+	// copy assignment
+	m = graph.m;
+    valid = true;
+}
+
+//template <typename Sample, typename Distance, typename WeightType, bool isDense, bool isSymmetric>
+//KNNGraph<Sample, Distance, WeightType, isDense, isSymmetric>::KNNGraph(Tree<Sample, Distance>& tree, size_t neighbors_num, size_t max_bruteforce_size, int max_iterations, double update_range)
+//    : Graph<WeightType, isDense, isSymmetric>(tree.size()), _neighbors_num(neighbors_num), _max_bruteforce_size(max_bruteforce_size), _max_iterations(max_iterations), _update_range(update_range)
+//{
+//	auto nodes = tree.get_all_nodes();
+//	std::vector<Sample> samples;
+//	std::vector<std::pair<size_t, size_t>> edgesPairs;
+//
+//	// iterate through the tree to collect edges	   	
+//	for (std::size_t i = 0; i < nodes.size(); ++i) 
+//	{
+//		samples.push_back(nodes[i]->data);
+//		auto children = nodes[i]->get_children();
+//		for (std::size_t j = 0; j < children.size(); ++j) 
+//		{
+//            edgesPairs.emplace_back(nodes[i]->get_ID(), children[j]->get_ID());
+//		}
+//	}
+//	
+//	_nodes = samples;
+//
+//	calculate_distance_matrix(samples);
+//
+//    buildEdges(edgesPairs);
+//
+//    valid = true;
+//}
+
+template <typename Sample, typename Distance, typename WeightType, bool isDense, bool isSymmetric>
+void KNNGraph<Sample, Distance, WeightType, isDense, isSymmetric>::calculate_distance_matrix(std::vector<Sample> samples)
+{
+	Distance distance;
+	
+	for (int i = 0; i < samples.size(); i++)
+	{
+		// take each node
+		auto i_point = samples[i];
+		std::vector<Distance::value_type> distances;
+		// then calculate distances for all other nodes
+		for (int j = 0; j < samples.size(); j++)
+		{
+			auto i_other_point = samples[j];
+			distances.push_back(distance(i_point, i_other_point));
+		}
+
+		_distance_matrix.push_back(distances);
+	}
+}
+
+template <typename Sample, typename Distance, typename WeightType, bool isDense, bool isSymmetric>
+void KNNGraph<Sample, Distance, WeightType, isDense, isSymmetric>::make_edge_pairs(std::vector<Sample> samples)
+{
+    m.resize(samples.size(), samples.size());	
+	
+	std::vector<int> ids(samples.size());
+    std::iota(ids.begin(), ids.end(), 0);
+
+	std::vector<std::pair<size_t, size_t>> edgesPairs;
+	double updated_percent = 1.0;
+	int iterations = 0;
+
+	// we iterate until we can found new edges (updated edges between nodes)
+	while(updated_percent > _update_range)
+	{
+		// create or update approximated knn graph
+		auto newEdgesPairs = random_pair_division(samples, ids, _max_bruteforce_size);
+					
+		// then update edge pair and check how many was updated
+		int was_size = edgesPairs.size();
+		for (int j = 0; j < newEdgesPairs.size(); j++)
+		{
+			bool already_exist = false;
+
+			for (int k = 0; k < edgesPairs.size(); k++)
+			{
+				if (edgesPairs[k] == newEdgesPairs[j] || (edgesPairs[k].first == newEdgesPairs[j].second && edgesPairs[k].second == newEdgesPairs[j].first))
+				{
+					already_exist = true;
+					break;
+				}
+			}
+
+			if (!already_exist)
+			{
+				edgesPairs.push_back(newEdgesPairs[j]);
+			}
+		}
+		if (was_size > 0)
+		{
+			updated_percent = (double)(edgesPairs.size() - was_size) / was_size;
+		}
+
+		iterations++;
+
+		if (iterations >= _max_iterations)
+		{
+			break;
+		}
+	}
+
+	// finish graph
+    buildEdges(edgesPairs);
+}
+
+template <typename Sample, typename Distance, typename WeightType, bool isDense, bool isSymmetric>
+void KNNGraph<Sample, Distance, WeightType, isDense, isSymmetric>::construct(std::vector<Sample> samples)
+{
+	calculate_distance_matrix(samples);
+
+    make_edge_pairs(samples);
+
+    valid = true;
+}
+
+
+template <typename Sample, typename Distance, typename WeightType, bool isDense, bool isSymmetric>
+std::vector<std::pair<size_t, size_t>> KNNGraph<Sample, Distance, WeightType, isDense, isSymmetric>::random_pair_division(std::vector<Sample> samples, std::vector<int> ids, int max_size)
+{
+	Distance d;
+	std::vector<Sample> A;
+	std::vector<Sample> B;
+	std::vector<int> A_ids;
+	std::vector<int> B_ids;
+    std::vector<std::pair<size_t, size_t>> edgesPairs;
+    std::vector<std::pair<size_t, size_t>> edgesPairsResult;
+
+	auto n = samples.size();
+
+	if (n > 0)
+	{
+		std::random_device rnd;
+		std::mt19937 mt(rnd());
+		std::uniform_int_distribution<int> dist(0, n - 1);
+
+		if (n <= max_size)
+		{
+			// conquer stage
+			edgesPairs = brute_force(samples, ids);
+		}
+		else
+		{
+			// divide stage
+
+			// take random nodes(samples)
+			auto a = samples[dist(mt)];
+			auto b = samples[dist(mt)];
+			Sample x;
+			// and divide all nodes to two groups, where each node is close to one of two initial points
+			for (int i = 0; i < n; i++)
+			{
+				x = samples[i];
+				if (d(x, a) < d(x, b))
+				{
+					A.push_back(x);
+					A_ids.push_back(ids[i]);
+				}
+				else
+				{
+					B.push_back(x);
+					B_ids.push_back(ids[i]);
+				}
+			}
+
+			// and recursively divide both groups again
+			edgesPairsResult = random_pair_division(A, A_ids, max_size);
+			edgesPairs.insert( edgesPairs.end(), edgesPairsResult.begin(), edgesPairsResult.end() );
+			edgesPairsResult = random_pair_division(B, B_ids, max_size);
+			edgesPairs.insert( edgesPairs.end(), edgesPairsResult.begin(), edgesPairsResult.end() );
+		}
+	}
+
+	return edgesPairs;
+}
+
+template <typename Sample, typename Distance, typename WeightType, bool isDense, bool isSymmetric>
+std::vector<std::pair<size_t, size_t>> KNNGraph<Sample, Distance, WeightType, isDense, isSymmetric>::brute_force(std::vector<Sample> samples, std::vector<int> ids)
+{	
+    std::vector<std::pair<size_t, size_t>> edgesPairs;
+
+	int update_count = 0;
+	
+	std::vector<std::vector<Distance::value_type>> distances;
+	std::vector<Distance::value_type> distance_row;
+	for (int i = 0; i < ids.size(); i++)
+	{
+		distance_row.clear();
+		for (int j = 0; j < ids.size(); j++)
+		{
+			distance_row.push_back(_distance_matrix[i][j]);
+		}
+		distances.push_back(distance_row);
+	}
+	
+    for (int i = 0; i < ids.size(); i++) 
+	{
+		auto idxs = sort_indexes(distances[i]);
+
+		for (int j = 0; j < idxs.size(); j++)
+		{
+			// omit first item because it is pair between the same item
+			if (j == 0)
+			{
+				continue;
+			}
+
+			// and break if we get '_neighbors_num' edges
+			if (j >= _neighbors_num)
+			{
+				break;
+			}
+			
+			bool already_exist = false;
+			
+			std::vector<int>::iterator max_index = std::max_element(ids.begin(), ids.end());
+			// here we keep number of edges from each node
+			std::vector<int> num_edges_by_node(ids[std::distance(ids.begin(), max_index)] + 1, 0);
+			for (int k = 0; k < edgesPairs.size(); k++)
+			{
+				num_edges_by_node[edgesPairs[k].first]++;
+				num_edges_by_node[edgesPairs[k].second]++;
+				if (edgesPairs[k] == std::pair<size_t, size_t>(ids[i], ids[idxs[j]]) || edgesPairs[k] == std::pair<size_t, size_t>(ids[idxs[j]], ids[i]))
+				{
+					already_exist = true;
+					break;
+				}
+				// if we want to keep neighbours strickt not more then _neighbors_num
+				if (_not_more_neighbors)
+				{
+					if (num_edges_by_node[edgesPairs[k].first] >= _neighbors_num || num_edges_by_node[edgesPairs[k].second] >= _neighbors_num )
+					{
+						already_exist = true;
+						break;
+					}
+				}
+			}
+
+			// add current node and closest to result
+			if (!already_exist)
+			{
+				edgesPairs.emplace_back(ids[i], ids[idxs[j]]);
+			}
+		}
+    }
+
+	return edgesPairs;
+}
+
+template <typename Sample, typename Distance, typename WeightType, bool isDense, bool isSymmetric>
+std::vector<int> KNNGraph<Sample, Distance, WeightType, isDense, isSymmetric>::gnnn_search(Sample query, int max_closest_num, int iterations, int num_greedy_moves, int num_expansions)
+{
+	std::vector<int> result;
+
+	// variables for choosen nodes during search
+	std::vector<Distance::value_type> choosen_distances;
+	std::vector<int> choosen_nodes;
+
+	// temp variables
+	std::vector<Distance::value_type> distances;
+	Distance::value_type distance;
+			
+	Distance distancer;
+	
+	// num_expansions should be less then k(neighbors_num) of the graph
+	if (num_expansions > _neighbors_num)
+	{
+		num_expansions = _neighbors_num;
+	}
+
+	// if params missed
+	if (num_expansions < 0)
+	{
+		num_expansions = _neighbors_num;
+	}
+	// if still negative
+	if (num_expansions < 0)
+	{
+		num_expansions = 1;
+	}
+	if (num_greedy_moves < 0)
+	{
+		// if param is missed we choose 20% of all nodes as number of mooves
+		num_greedy_moves = round(_nodes.size() * 0.2);
+	}
+
+	//
+
+	if (_nodes.size() > 0)
+	{
+		for (int i = 0; i < iterations; i++) 
+		{
+			std::random_device rnd;
+			std::mt19937 mt(rnd());
+			std::uniform_int_distribution<int> dist(0, _nodes.size() - 1);
+		
+			// get initial random node from the graph
+			int checking_node = dist(mt);
+			int prev_node = -1;
+			int new_node;
+		
+			// walk from initial node on distance 'num_greedy_moves' steps
+			for (int j = 0; j < num_greedy_moves; j++) 
+			{
+				distances.clear();
+				// 0 index is for node itself, 1 - is first circle of neighbours
+				auto neighbours = getNeighbours(checking_node, 1)[1];
+
+				// get first num_expansions neighbours for the checking node and calculate distances to the query
+				for (int p = 0; p < num_expansions; p++) 
+				{
+					if (p < neighbours.size())
+					{
+						distance = distancer(_nodes[neighbours[p]], query);
+						distances.push_back(distance);
+					
+						if (std::find(choosen_nodes.begin(), choosen_nodes.end(), neighbours[p]) == choosen_nodes.end())
+						{
+							choosen_distances.push_back(distance);
+							choosen_nodes.push_back(neighbours[p]);
+						}
+					}
+				}
+			
+				std::vector<Distance::value_type>::iterator min_index = std::min_element(distances.begin(), distances.end());
+				new_node = neighbours[std::distance(distances.begin(), min_index)];
+				// if we back to the visited node then we fall in loop and search is complete
+				if (new_node == prev_node)
+				{
+					break;
+				}
+				prev_node = checking_node;
+				checking_node = new_node;
+
+			
+			}
+		}
+			
+		// sort distances and return corresopnding nodes from choosen
+		auto idxs = sort_indexes(choosen_distances);
+		for (int i = 0; i < max_closest_num; i++)
+		{
+			if (i < idxs.size())
+			{
+				result.push_back(choosen_nodes[idxs[i]]);
+			}
+		}
+
+	}
+
+	return result;
+}
+
+
+template <typename Sample, typename Distance, typename WeightType, bool isDense, bool isSymmetric>
+template <typename T1>
+std::vector<size_t> KNNGraph<Sample, Distance, WeightType, isDense, isSymmetric>::sort_indexes(const std::vector<T1> &v) 
+{
+	// initialize original index locations
+	std::vector<size_t> idx(v.size());
+	std::iota(idx.begin(), idx.end(), 0);
+
+	// sort indexes based on comparing values in v
+	std::sort(idx.begin(), idx.end(),
+		[&v](size_t i1, size_t i2) {return v[i1] < v[i2];});
+
+	return idx;
+}
+
 // Graph factory based on Blaze matrix of 4 allowed types
 
 template <class ValueType>
diff --git a/tests/mapping_tests/CMakeLists.txt b/tests/mapping_tests/CMakeLists.txt
index 95973aff..bd8d2c46 100644
--- a/tests/mapping_tests/CMakeLists.txt
+++ b/tests/mapping_tests/CMakeLists.txt
@@ -2,7 +2,8 @@ cmake_minimum_required(VERSION 3.10)
 
 enable_testing()
 
-set(Boost_DEBUG on)
+set(Boost_DEBUG on)  
+find_package(Boost)
 find_package(LAPACK)
 
 project(mapping_tests)
@@ -16,13 +17,11 @@ endif(MSVC)
 
 set(BOOST_LINK_DIRS "${BOOST_ROOT}/stage/lib")
 
-# mapping tests
-
-add_executable(mapping_tests mapping_tests.cpp)
+# mapping testsadd_executable(mapping_tests mapping_tests.cpp)
 
 target_include_directories(mapping_tests PRIVATE ${Boost_INCLUDE_DIRS})
+target_link_libraries(mapping_tests ${Boost_LIBRARIES})
 target_link_directories(mapping_tests PRIVATE ${BOOST_LINK_DIRS})
-target_link_libraries(mapping_tests PRIVATE ${Boost_LIBRARIES})
 
 if(UNIX)
     set(THREADS_PREFER_PTHREAD_FLAG ON)
diff --git a/tests/mapping_tests/mapping_tests.cpp b/tests/mapping_tests/mapping_tests.cpp
index f019534a..8ee3da5e 100644
--- a/tests/mapping_tests/mapping_tests.cpp
+++ b/tests/mapping_tests/mapping_tests.cpp
@@ -12,7 +12,7 @@
 #include "modules/utils/graph.hpp"
 
 #define BOOST_TEST_MODULE Main
-#define BOOST_TEST_DYN_LINK
+//#define BOOST_TEST_DYN_LINK
 
 #include <boost/test/unit_test.hpp>
 
diff --git a/tests/space_tests/CMakeLists.txt b/tests/space_tests/CMakeLists.txt
index 7b45704f..8b71bd27 100644
--- a/tests/space_tests/CMakeLists.txt
+++ b/tests/space_tests/CMakeLists.txt
@@ -1,12 +1,24 @@
+cmake_minimum_required(VERSION 3.10)
+
+enable_testing()
+
 set(BOOST_TEST_LIB ${Boost_LIBRARIES})
-message(${BOOST_TEST_LIB})
+# message(${BOOST_TEST_LIB})
 find_package (Boost COMPONENTS serialization system REQUIRED)
+find_package(LAPACK)
+
 include_directories(
     ${PROJECT_SOURCE_DIR}
+	${PROJECT_SOURCE_DIR}/../..
     ${Boost_INCLUDE_DIRS}
     )
 
-
+# for Visual Studio
+if(MSVC) 
+	set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} /std:c++17")
+endif(MSVC)
+	
+set(BOOST_LINK_DIRS "${BOOST_ROOT}/stage/lib")
 
 if(UNIX)
     set(THREADS_PREFER_PTHREAD_FLAG ON)
@@ -25,15 +37,30 @@ foreach(testSrc ${TEST_SRCS})
 
         #link to Boost libraries AND your targets and dependencies
 	if(WIN32)
-		target_compile_options(${testName} PRIVATE /D BOOST_ALL_NO_LIB=1)
+		#target_compile_options(${testName} PRIVATE /D BOOST_ALL_NO_LIB=1)
 #		target_compile_options(${testName} PRIVATE /DBOOST_USE_WINDOWS_H)
 	endif()
-        target_link_libraries(${testName} ${Boost_LIBRARIES})
-        target_link_libraries(${testName} ${BOOST_TEST_LIB})
+		target_include_directories(${testName} PRIVATE ${Boost_INCLUDE_DIRS})
+		target_link_directories(${testName} PRIVATE ${BOOST_LINK_DIRS})
+        target_link_libraries(${testName} PRIVATE ${Boost_LIBRARIES})
+        target_link_libraries(${testName} PRIVATE ${BOOST_TEST_LIB})
+
+	if(LAPACK_LIBRARIES)
+		target_link_libraries(${testName} PRIVATE ${LAPACK_LIBRARIES})
+	endif(LAPACK_LIBRARIES)
+ 
+	if(CMAKE_SYSTEM_NAME MATCHES Windows)
+		target_link_directories(${testName} PRIVATE ${PROJECT_SOURCE_DIR}/../../mkl/lib/intel64_win)
+		target_link_libraries(${testName} PRIVATE mkl_core.lib mkl_sequential.lib mkl_intel_lp64.lib)
+		target_link_libraries(${testName} PRIVATE ${Boost_LIBRARIES})
+	endif()
+
 	if(UNIX)
 	    target_link_libraries(${testName} Threads::Threads)
+	set_target_properties(mapping_tests PROPERTIES CXX_STANDARD 17)
 	endif(UNIX)
-        #I like to move testing binaries into a testBin directory
+
+        #I like to move testing binaries into a testBin directory 
 #        set_target_properties(${testName} PROPERTIES RUNTIME_OUTPUT_DIRECTORY  ${CMAKE_CURRENT_SOURCE_DIR})
         #Finally add it to test execution - 
         #Notice the WORKING_DIRECTORY and COMMAND
diff --git a/tests/space_tests/README.md b/tests/space_tests/README.md
new file mode 100644
index 00000000..ca20eb6f
--- /dev/null
+++ b/tests/space_tests/README.md
@@ -0,0 +1,33 @@
+
+---
+
+## Run
+*You need STL and C++17 support to compile.*
+
+#### Using CMake
+
+_Windows_
+
+- First of all you need to install [Intel MKL](https://software.intel.com/en-us/mkl)
+- Then copy `mkl` (usually from `C:\Program Files (x86)\IntelSWTools\compilers_and_libraries\windows`) folder to the root of the project. 
+- Then you can run cmake to create project (links to `mkl` already inside `CMakeLists.txt`)
+```bash
+mkdir build
+cd build
+cmake .. -A x64 -DBoost_NAMESPACE="libboost" -DBoost_COMPILER="-vc141" -DBOOST_ROOT=C:/boost_1_70_0/
+```
+
+Where `DBoost_NAMESPACE` and `DBoost_COMPILER` params point to suffixes and postfixes of the search libs in the `stage/lib` folder of the boost root directory. 
+And `DBOOST_ROOT` point to the boost root folder. 
+
+Then open solution in the Microsoft Visual Studio.
+
+_Linux_
+
+Just run cmake
+```bash
+mkdir build
+cd build
+cmake ..
+make
+```
\ No newline at end of file
diff --git a/tests/space_tests/knn_graph_test.cpp b/tests/space_tests/knn_graph_test.cpp
new file mode 100644
index 00000000..80c0aafd
--- /dev/null
+++ b/tests/space_tests/knn_graph_test.cpp
@@ -0,0 +1,313 @@
+/*
+This Source Code Form is subject to the terms of the Mozilla Public
+License, v. 2.0. If a copy of the MPL was not distributed with this
+file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+Copyright (c) 2019 Panda Team
+*/
+
+#define BOOST_TEST_MODULE knn_graph_test
+#include <boost/test/unit_test.hpp>
+
+#include <iostream>
+#include "modules/utils/graph.hpp"
+#include "modules/distance.hpp"
+
+BOOST_AUTO_TEST_CASE(valid_params)
+{	
+    size_t neighbors_num = 3;
+    size_t search_neighbors_num = 3;
+		
+    std::vector<std::vector<double>> table = {
+        { 0, 1 },
+        { 1, 1 },
+        { 2, 2 },
+        { 3, 3 },
+        { 4, 3 },
+        { 5, 3 },
+        { 4, 6 },
+        { 5, 1 },
+        { 4, 1 },
+        { 3, 2 },
+        { 0, 3 },
+        { 1, 3 },
+        { 2, 3 },
+        { 6, 6 },
+        { 7, 6 },
+        { 0, 2 },
+        { 0, 9 },
+        { 0, 4 },
+        { 0, 5 },
+        { 0, 6 },
+        { 0, 7 },
+        { 0, 8 },
+    };
+	
+	std::vector<double> query = { 7, 7 };
+	
+	///	
+
+	// valid params
+    auto g_0 = metric::KNNGraph<std::vector<double>, metric::Euclidian<double>>(table, neighbors_num, 2.5 * neighbors_num);
+
+	// nothing bad should be happened
+    std::cout << std::endl;
+    std::cout << " --- valid params --- " << std::endl;
+    std::cout << "gnn search:" << std::endl;
+	auto found_0 = g_0.gnnn_search(query, search_neighbors_num);
+    std::sort(found_0.begin(), found_0.end());
+	print_vector(found_0);
+	BOOST_CHECK(found_0.size() == 3);
+	// we expect only 4, 5, 6, 13 or 14 nodes (depends on random)
+	BOOST_CHECK(found_0[0] == 4 || found_0[0] == 5 || found_0[0] == 6 || found_0[0] == 13 || found_0[0] == 14 );
+	BOOST_CHECK(found_0[1] == 4 || found_0[1] == 5 || found_0[1] == 6 || found_0[1] == 13 || found_0[1] == 14 );
+	BOOST_CHECK(found_0[2] == 4 || found_0[2] == 5 || found_0[2] == 6 || found_0[2] == 13 || found_0[2] == 14 );
+    std::cout << std::endl;
+}
+
+BOOST_AUTO_TEST_CASE(negative_params)
+{
+    size_t neighbors_num = 3;
+    size_t search_neighbors_num = 3;
+		
+    std::vector<std::vector<double>> table = {
+        { 0, 1 },
+        { 1, 1 },
+        { 2, 2 },
+        { 3, 3 },
+        { 4, 3 },
+        { 5, 3 },
+        { 4, 6 },
+        { 5, 1 },
+        { 4, 1 },
+        { 3, 2 },
+        { 0, 3 },
+        { 1, 3 },
+        { 2, 3 },
+        { 6, 6 },
+        { 7, 6 },
+        { 0, 2 },
+        { 0, 9 },
+        { 0, 4 },
+        { 0, 5 },
+        { 0, 6 },
+        { 0, 7 },
+        { 0, 8 },
+    };
+	
+	std::vector<double> query = { 7, 7 };
+
+	///	
+
+	// negative params
+    auto g_1 = metric::KNNGraph<std::vector<double>, metric::Euclidian<double>>(table, -neighbors_num, -2.5 * neighbors_num);
+
+	// nothing bad should be happened
+    std::cout << std::endl;
+    std::cout << " --- negative params --- " << std::endl;
+    std::cout << "gnn search:" << std::endl;
+	auto found_1 = g_1.gnnn_search(query, search_neighbors_num);
+	print_vector(found_1);
+	BOOST_CHECK(found_1.size() > 0);
+    std::cout << std::endl;
+
+	// negative params for search, nothing bad should be happened too
+    std::cout << std::endl;
+    std::cout << "gnn search:" << std::endl;
+	found_1 = g_1.gnnn_search(query, -search_neighbors_num, -10);
+	print_vector(found_1);
+	BOOST_CHECK(found_1.size() == 0);
+    std::cout << std::endl;
+	
+}
+
+BOOST_AUTO_TEST_CASE(params_more_than_nodes_exists)
+{
+    size_t neighbors_num = 3;
+    size_t search_neighbors_num = 3;
+		
+    std::vector<std::vector<double>> table = {
+        { 0, 1 },
+        { 1, 1 },
+        { 2, 2 },
+        { 3, 3 },
+        { 4, 3 },
+        { 5, 3 },
+        { 4, 6 },
+        { 5, 1 },
+        { 4, 1 },
+        { 3, 2 },
+        { 0, 3 },
+        { 1, 3 },
+        { 2, 3 },
+        { 6, 6 },
+        { 7, 6 },
+        { 0, 2 },
+        { 0, 9 },
+        { 0, 4 },
+        { 0, 5 },
+        { 0, 6 },
+        { 0, 7 },
+        { 0, 8 },
+    };
+	
+	std::vector<double> query = { 7, 7 };
+
+	///
+
+	// params more than nodes exists
+    auto g_2 = metric::KNNGraph<std::vector<double>, metric::Euclidian<double>>(table, table.size() * 10, table.size() * 25);
+
+	// nothing bad should be happened
+    std::cout << std::endl;
+    std::cout << " --- params more than nodes exists --- " << std::endl;
+    std::cout << "gnn search:" << std::endl;
+	auto found_2 = g_2.gnnn_search(query, table.size() * 50);
+	print_vector(found_2);
+	BOOST_CHECK(found_2.size() > 0);
+    std::cout << std::endl;
+}
+
+BOOST_AUTO_TEST_CASE(empty_dataset)
+{
+    size_t neighbors_num = 3;
+    size_t search_neighbors_num = 3;
+		
+    std::vector<std::vector<double>> table = {
+        { 0, 1 },
+        { 1, 1 },
+        { 2, 2 },
+        { 3, 3 },
+        { 4, 3 },
+        { 5, 3 },
+        { 4, 6 },
+        { 5, 1 },
+        { 4, 1 },
+        { 3, 2 },
+        { 0, 3 },
+        { 1, 3 },
+        { 2, 3 },
+        { 6, 6 },
+        { 7, 6 },
+        { 0, 2 },
+        { 0, 9 },
+        { 0, 4 },
+        { 0, 5 },
+        { 0, 6 },
+        { 0, 7 },
+        { 0, 8 },
+    };
+	
+	std::vector<double> query = { 7, 7 };
+
+	///
+
+	// empty dataset
+	auto g_3 = metric::KNNGraph<std::vector<double>, metric::Euclidian<double>>({}, neighbors_num, 2.5 * neighbors_num);
+
+	// nothing bad should be happened
+    std::cout << std::endl;
+    std::cout << " --- empty dataset --- " << std::endl;
+    std::cout << "gnn search:" << std::endl;
+	auto found_3 = g_3.gnnn_search(query, search_neighbors_num);
+	print_vector(found_3);
+	BOOST_CHECK(found_3.size() == 0);
+    std::cout << std::endl;
+}
+
+BOOST_AUTO_TEST_CASE(dataset_and_query_with_different_dimensions)
+{
+    size_t neighbors_num = 3;
+    size_t search_neighbors_num = 3;
+		
+    std::vector<std::vector<double>> table = {
+        { 0, 1 },
+        { 1, 1 },
+        { 2, 2 },
+        { 3, 3 },
+        { 4, 3 },
+        { 5, 3 },
+        { 4, 6 },
+        { 5, 1 },
+        { 4, 1 },
+        { 3, 2 },
+        { 0, 3 },
+        { 1, 3 },
+        { 2, 3 },
+        { 6, 6 },
+        { 7, 6 },
+        { 0, 2 },
+        { 0, 9 },
+        { 0, 4 },
+        { 0, 5 },
+        { 0, 6 },
+        { 0, 7 },
+        { 0, 8 },
+    };
+	
+	std::vector<double> query = { 7, 7 };
+
+	///
+
+	// dataset and query with different dimensions
+	auto g_4 = metric::KNNGraph<std::vector<double>, metric::Euclidian<double>>({}, neighbors_num, 2.5 * neighbors_num);
+
+	// nothing bad should be happened
+    std::cout << std::endl;
+    std::cout << " --- dataset and query with different dimensions --- " << std::endl;
+    std::cout << "gnn search:" << std::endl;
+	auto found_4 = g_4.gnnn_search({ 7, 5, 7, 5, 7, 5, 7, 5 }, search_neighbors_num);
+	print_vector(found_4);
+	BOOST_CHECK(found_4.size() == 0);
+    std::cout << std::endl;
+}
+
+BOOST_AUTO_TEST_CASE(copy_constructor)
+{
+    size_t neighbors_num = 3;
+    size_t search_neighbors_num = 3;
+		
+    std::vector<std::vector<double>> table = {
+        { 0, 1 },
+        { 1, 1 },
+        { 2, 2 },
+        { 3, 3 },
+        { 4, 3 },
+        { 5, 3 },
+        { 4, 6 },
+        { 5, 1 },
+        { 4, 1 },
+        { 3, 2 },
+        { 0, 3 },
+        { 1, 3 },
+        { 2, 3 },
+        { 6, 6 },
+        { 7, 6 },
+        { 0, 2 },
+        { 0, 9 },
+        { 0, 4 },
+        { 0, 5 },
+        { 0, 6 },
+        { 0, 7 },
+        { 0, 8 },
+    };
+	
+	std::vector<double> query = { 7, 7 };
+
+	///
+	
+	//	
+    auto g_0 = metric::KNNGraph<std::vector<double>, metric::Euclidian<double>>(table, neighbors_num, 2.5 * neighbors_num);
+	auto found_0 = g_0.gnnn_search(query, search_neighbors_num);
+
+    auto g_5 = metric::KNNGraph(g_0);
+
+	// nothing bad should be happened
+    std::cout << std::endl;
+    std::cout << " --- nothing bad should be happened --- " << std::endl;
+    std::cout << "gnn search:" << std::endl;
+	auto found_5 = g_5.gnnn_search(query, search_neighbors_num);
+	print_vector(found_5);
+	assert(found_5.size() == found_0.size());
+    std::cout << std::endl;
+}