added possibility of negative values in persistance diagram

1 files changed, 41 insertions, 16 deletions

diff --git a/src/Bottleneck_distance/include/gudhi/Persistence_graph.h b/src/Bottleneck_distance/include/gudhi/Persistence_graph.h
index 33f03b9c..9b663cc2 100644
--- a/src/Bottleneck_distance/include/gudhi/Persistence_graph.h
+++ b/src/Bottleneck_distance/include/gudhi/Persistence_graph.h
@@ -20,6 +20,7 @@
 #include <vector>
 #include <algorithm>
 #include <limits>  // for numeric_limits
+#include <cmath>
 
 namespace Gudhi {
 
@@ -70,27 +71,51 @@ Persistence_graph::Persistence_graph(const Persistence_diagram1 &diag1,
     : u(), v(), b_alive(0.) {
   std::vector<double> u_alive;
   std::vector<double> v_alive;
+  std::vector<double> u_nalive;
+  std::vector<double> v_nalive;
+  int u_inf = 0;
+  int v_inf = 0;
+  double inf = std::numeric_limits<double>::infinity();
+  double neginf = -1 * inf;
+
   for (auto it = std::begin(diag1); it != std::end(diag1); ++it) {
-    if (std::get<1>(*it) == std::numeric_limits<double>::infinity())
-      u_alive.push_back(std::get<0>(*it));
-    else if (std::get<1>(*it) - std::get<0>(*it) > e)
-      u.push_back(Internal_point(std::get<0>(*it), std::get<1>(*it), u.size()));
+	if (std::get<0>(*it) != inf && std::get<1>(*it) != neginf){
+	  if (std::get<0>(*it) == neginf && std::get<1>(*it) == inf)
+		u_inf++;
+	  else if (std::get<0>(*it) == neginf)
+		u_nalive.push_back(std::get<1>(*it));
+	  else if (std::get<1>(*it) == inf)
+		u_alive.push_back(std::get<0>(*it));
+	  else if (std::get<1>(*it) - std::get<0>(*it) > e)
+		u.push_back(Internal_point(std::get<0>(*it), std::get<1>(*it), u.size()));
+	}
   }
   for (auto it = std::begin(diag2); it != std::end(diag2); ++it) {
-    if (std::get<1>(*it) == std::numeric_limits<double>::infinity())
-      v_alive.push_back(std::get<0>(*it));
-    else if (std::get<1>(*it) - std::get<0>(*it) > e)
-      v.push_back(Internal_point(std::get<0>(*it), std::get<1>(*it), v.size()));
+	if (std::get<0>(*it) != inf && std::get<1>(*it) != neginf){
+	  if (std::get<0>(*it) == neginf && std::get<1>(*it) == inf)
+		v_inf++;
+	  else if (std::get<0>(*it) == neginf)
+		v_nalive.push_back(std::get<1>(*it));
+	  else if (std::get<1>(*it) == inf)
+		v_alive.push_back(std::get<0>(*it));
+	  else if (std::get<1>(*it) - std::get<0>(*it) > e)
+		v.push_back(Internal_point(std::get<0>(*it), std::get<1>(*it), v.size()));
+	}
   }
   if (u.size() < v.size())
     swap(u, v);
-  std::sort(u_alive.begin(), u_alive.end());
-  std::sort(v_alive.begin(), v_alive.end());
-  if (u_alive.size() != v_alive.size()) {
+
+  if (u_alive.size() != v_alive.size() || u_nalive.size() != v_nalive.size() || u_inf != v_inf) {
     b_alive = std::numeric_limits<double>::infinity();
   } else {
+	std::sort(u_alive.begin(), u_alive.end());
+	std::sort(v_alive.begin(), v_alive.end());
+	std::sort(u_nalive.begin(), u_nalive.end());
+	std::sort(v_nalive.begin(), v_nalive.end());
     for (auto it_u = u_alive.cbegin(), it_v = v_alive.cbegin(); it_u != u_alive.cend(); ++it_u, ++it_v)
       b_alive = (std::max)(b_alive, std::fabs(*it_u - *it_v));
+	for (auto it_u = u_nalive.cbegin(), it_v = v_nalive.cbegin(); it_u != u_nalive.cend(); ++it_u, ++it_v)
+	  b_alive = (std::max)(b_alive, std::fabs(*it_u - *it_v));
   }
 }
 
@@ -114,7 +139,7 @@ inline int Persistence_graph::corresponding_point_in_v(int u_point_index) const
 
 inline double Persistence_graph::distance(int u_point_index, int v_point_index) const {
   if (on_the_u_diagonal(u_point_index) && on_the_v_diagonal(v_point_index))
-    return 0.;
+	return 0.;
   Internal_point p_u = get_u_point(u_point_index);
   Internal_point p_v = get_v_point(v_point_index);
   return (std::max)(std::fabs(p_u.x() - p_v.x()), std::fabs(p_u.y() - p_v.y()));
@@ -132,9 +157,9 @@ inline std::vector<double> Persistence_graph::sorted_distances() const {
   std::vector<double> distances;
   distances.push_back(0.);  // for empty diagrams
   for (int u_point_index = 0; u_point_index < size(); ++u_point_index) {
-    distances.push_back(distance(u_point_index, corresponding_point_in_v(u_point_index)));
-    for (int v_point_index = 0; v_point_index < size(); ++v_point_index)
-      distances.push_back(distance(u_point_index, v_point_index));
+	distances.push_back(distance(u_point_index, corresponding_point_in_v(u_point_index)));
+	for (int v_point_index = 0; v_point_index < size(); ++v_point_index)
+	  distances.push_back(distance(u_point_index, v_point_index));
   }
 #ifdef GUDHI_USE_TBB
   tbb::parallel_sort(distances.begin(), distances.end());
@@ -146,7 +171,7 @@ inline std::vector<double> Persistence_graph::sorted_distances() const {
 
 inline Internal_point Persistence_graph::get_u_point(int u_point_index) const {
   if (!on_the_u_diagonal(u_point_index))
-    return u.at(u_point_index);
+	return u.at(u_point_index);
   Internal_point projector = v.at(corresponding_point_in_v(u_point_index));
   double m = (projector.x() + projector.y()) / 2.;
   return Internal_point(m, m, u_point_index);