alpha_complex_persistence utils can use weighted points. Update Alpha_kernel_d_unit_test to work

4 files changed, 84 insertions, 54 deletions

diff --git a/src/Alpha_complex/test/Alpha_kernel_d_unit_test.cpp b/src/Alpha_complex/test/Alpha_kernel_d_unit_test.cpp
index 70b65aad..192834b3 100644
--- a/src/Alpha_complex/test/Alpha_kernel_d_unit_test.cpp
+++ b/src/Alpha_complex/test/Alpha_kernel_d_unit_test.cpp
@@ -102,31 +102,3 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_kernel_d_sphere, TestedKernel, list_of_kerne
   std::clog << "Squared radius is " << w_sq_rd << std::endl;
   GUDHI_TEST_FLOAT_EQUALITY_CHECK(cast_to_double(w_sphere.weight()), cast_to_double(w_sq_rd));
 }
-
-
-BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_kernel_d_distance, TestedKernel, list_of_kernel_variants) {
-  using Unweighted_kernel = Gudhi::alpha_complex::Alpha_kernel_d<TestedKernel, false>;
-
-  std::vector<double> p0 {1., 0., 0., 0.};
-  std::vector<double> p1 {0., 1., 0., 0.};
-
-  using Point_d = typename Unweighted_kernel::Point_d;
-  Unweighted_kernel kernel;
-  auto dist_01 = kernel.get_squared_distance(Point_d(p0.begin(), p0.end()), Point_d(p1.begin(), p1.end()));
-  std::clog << "Distance is " << dist_01 << std::endl;
-
-  using Weighted_kernel = Gudhi::alpha_complex::Alpha_kernel_d<TestedKernel, true>;
-
-  using Weighted_point_d = typename Weighted_kernel::Weighted_point_d;
-  using Bare_point_d = typename Weighted_kernel::Bare_point_d;
-  std::vector<Weighted_point_d> w_pts;
-
-  Weighted_kernel w_kernel;
-  auto w_dist_01 = w_kernel.get_squared_distance(Weighted_point_d(Bare_point_d(p0.begin(), p0.end()), 0.),
-                                                 Weighted_point_d(Bare_point_d(p1.begin(), p1.end()), 0.));
-  std::clog << "Distance is " << w_dist_01 << std::endl;
-
-  CGAL::NT_converter<typename Weighted_kernel::FT, double> cast_to_double;
-  // The results shall be the same with weights = 0.
-  GUDHI_TEST_FLOAT_EQUALITY_CHECK(cast_to_double(dist_01), cast_to_double(w_dist_01));
-}
diff --git a/src/Alpha_complex/test/Weighted_alpha_complex_unit_test.cpp b/src/Alpha_complex/test/Weighted_alpha_complex_unit_test.cpp
index b4fc76de..57a57058 100644
--- a/src/Alpha_complex/test/Weighted_alpha_complex_unit_test.cpp
+++ b/src/Alpha_complex/test/Weighted_alpha_complex_unit_test.cpp
@@ -21,6 +21,7 @@
 #include <random>
 
 #include <gudhi/Alpha_complex.h>
+#include <gudhi/Alpha_complex_3d.h>
 #include <gudhi/Simplex_tree.h>
 #include <gudhi/Unitary_tests_utils.h>
 
@@ -32,7 +33,6 @@ typedef CGAL::Epeck_d< CGAL::Dimension_tag<4> > Exact_kernel_s;
 typedef CGAL::Epick_d< CGAL::Dynamic_dimension_tag > Inexact_kernel_d;
 // Use static dimension_tag for the user not to be able to set dimension
 typedef CGAL::Epick_d< CGAL::Dimension_tag<4> > Inexact_kernel_s;
-// The triangulation uses the default instantiation of the TriangulationDataStructure template parameter
 
 typedef boost::mpl::list<Exact_kernel_d, Exact_kernel_s, Inexact_kernel_d, Inexact_kernel_s> list_of_kernel_variants;
 
diff --git a/src/Alpha_complex/utilities/alpha_complex_persistence.cpp b/src/Alpha_complex/utilities/alpha_complex_persistence.cpp
index e17831d9..e86b34e2 100644
--- a/src/Alpha_complex/utilities/alpha_complex_persistence.cpp
+++ b/src/Alpha_complex/utilities/alpha_complex_persistence.cpp
@@ -17,19 +17,82 @@
 #include <gudhi/Persistent_cohomology.h>
 // to construct a simplex_tree from alpha complex
 #include <gudhi/Simplex_tree.h>
+#include <gudhi/Points_off_io.h>
 
 #include <iostream>
 #include <string>
 #include <limits>  // for numeric_limits
+#include <vector>
+#include <fstream>
 
 using Simplex_tree = Gudhi::Simplex_tree<>;
 using Filtration_value = Simplex_tree::Filtration_value;
 
 void program_options(int argc, char *argv[], std::string &off_file_points, bool &exact, bool &fast,
-                     std::string &output_file_diag, Filtration_value &alpha_square_max_value,
+                     std::string &weight_file, std::string &output_file_diag, Filtration_value &alpha_square_max_value,
                      int &coeff_field_characteristic, Filtration_value &min_persistence);
 
+template<class Point_d>
+std::vector<Point_d> read_off(const std::string &off_file_points) {
+  Gudhi::Points_off_reader<Point_d> off_reader(off_file_points);
+  if (!off_reader.is_valid()) {
+    std::cerr << "Alpha_complex - Unable to read file " << off_file_points << "\n";
+    exit(-1);  // ----- >>
+  }
+  return off_reader.get_point_cloud();
+}
+
+std::vector<double> read_weight_file(const std::string &weight_file) {
+  std::vector<double> weights;
+  // Read weights information from file
+  std::ifstream weights_ifstr(weight_file);
+  if (weights_ifstr.good()) {
+    double weight = 0.0;
+    // Attempt read the weight in a double format, return false if it fails
+    while (weights_ifstr >> weight) {
+      weights.push_back(weight);
+    }
+  } else {
+    std::cerr << "Unable to read weights file " << weight_file << std::endl;
+    exit(-1);
+  }
+  return weights;
+}
+
+template<class Kernel>
+Simplex_tree create_simplex_tree(const std::string &off_file_points, const std::string &weight_file,
+                                 bool exact_version, Filtration_value alpha_square_max_value) {
+  Simplex_tree stree;
+  auto points = read_off<typename Kernel::Point_d>(off_file_points);
+
+  if (weight_file != std::string()) {
+    std::vector<double> weights = read_weight_file(weight_file);
+    if (points.size() != weights.size()) {
+      std::cerr << "Alpha_complex - Inconsistency between number of points (" << points.size()
+                << ") and number of weights (" << weights.size() << ")" << "\n";
+      exit(-1);  // ----- >>
+    }
+    // Init of an alpha complex from an OFF file
+    Gudhi::alpha_complex::Alpha_complex<Kernel, true> alpha_complex_from_file(points, weights);
+
+    if (!alpha_complex_from_file.create_complex(stree, alpha_square_max_value, exact_version)) {
+      std::cerr << "Alpha complex simplicial complex creation failed." << std::endl;
+      exit(-1);
+    }
+  } else {
+    // Init of an alpha complex from an OFF file
+    Gudhi::alpha_complex::Alpha_complex<Kernel> alpha_complex_from_file(points);
+
+    if (!alpha_complex_from_file.create_complex(stree, alpha_square_max_value, exact_version)) {
+      std::cerr << "Alpha complex simplicial complex creation failed." << std::endl;
+      exit(-1);
+    }
+  }
+  return stree;
+}
+
 int main(int argc, char **argv) {
+  std::string weight_file;
   std::string off_file_points;
   std::string output_file_diag;
   bool exact_version = false;
@@ -38,48 +101,34 @@ int main(int argc, char **argv) {
   int coeff_field_characteristic;
   Filtration_value min_persistence;
 
-  program_options(argc, argv, off_file_points, exact_version, fast_version, output_file_diag, alpha_square_max_value,
-                  coeff_field_characteristic, min_persistence);
+  program_options(argc, argv, off_file_points, exact_version, fast_version, weight_file, output_file_diag,
+                  alpha_square_max_value, coeff_field_characteristic, min_persistence);
 
   if ((exact_version) && (fast_version)) {
     std::cerr << "You cannot set the exact and the fast version." << std::endl;
     exit(-1);
   }
 
-  Simplex_tree simplex;
+  Simplex_tree stree;
   if (fast_version) {
     // WARNING : CGAL::Epick_d is fast but not safe (unlike CGAL::Epeck_d)
     // (i.e. when the points are on a grid)
     using Fast_kernel = CGAL::Epick_d<CGAL::Dynamic_dimension_tag>;
-
-    // Init of an alpha complex from an OFF file
-    Gudhi::alpha_complex::Alpha_complex<Fast_kernel> alpha_complex_from_file(off_file_points);
-
-    if (!alpha_complex_from_file.create_complex(simplex, alpha_square_max_value)) {
-      std::cerr << "Fast Alpha complex simplicial complex creation failed." << std::endl;
-      exit(-1);
-    }
+    stree = create_simplex_tree<Fast_kernel>(off_file_points, weight_file, exact_version, alpha_square_max_value);
   } else {
     using Kernel = CGAL::Epeck_d<CGAL::Dynamic_dimension_tag>;
-
-    // Init of an alpha complex from an OFF file
-    Gudhi::alpha_complex::Alpha_complex<Kernel> alpha_complex_from_file(off_file_points);
-
-    if (!alpha_complex_from_file.create_complex(simplex, alpha_square_max_value, exact_version)) {
-      std::cerr << "Alpha complex simplicial complex creation failed." << std::endl;
-      exit(-1);
-    }
+    stree = create_simplex_tree<Kernel>(off_file_points, weight_file, exact_version, alpha_square_max_value);
   }
   // ----------------------------------------------------------------------------
   // Display information about the alpha complex
   // ----------------------------------------------------------------------------
-  std::clog << "Simplicial complex is of dimension " << simplex.dimension() << " - " << simplex.num_simplices()
-            << " simplices - " << simplex.num_vertices() << " vertices." << std::endl;
+  std::clog << "Simplicial complex is of dimension " << stree.dimension() << " - " << stree.num_simplices()
+            << " simplices - " << stree.num_vertices() << " vertices." << std::endl;
 
-  std::clog << "Simplex_tree dim: " << simplex.dimension() << std::endl;
+  std::clog << "Simplex_tree dim: " << stree.dimension() << std::endl;
   // Compute the persistence diagram of the complex
   Gudhi::persistent_cohomology::Persistent_cohomology<Simplex_tree, Gudhi::persistent_cohomology::Field_Zp> pcoh(
-      simplex);
+      stree);
   // initializes the coefficient field for homology
   pcoh.init_coefficients(coeff_field_characteristic);
 
@@ -98,7 +147,7 @@ int main(int argc, char **argv) {
 }
 
 void program_options(int argc, char *argv[], std::string &off_file_points, bool &exact, bool &fast,
-                     std::string &output_file_diag, Filtration_value &alpha_square_max_value,
+                     std::string &weight_file, std::string &output_file_diag, Filtration_value &alpha_square_max_value,
                      int &coeff_field_characteristic, Filtration_value &min_persistence) {
   namespace po = boost::program_options;
   po::options_description hidden("Hidden options");
@@ -111,6 +160,8 @@ void program_options(int argc, char *argv[], std::string &off_file_points, bool
       "To activate exact version of Alpha complex (default is false, not available if fast is set)")(
       "fast,f", po::bool_switch(&fast),
       "To activate fast version of Alpha complex (default is false, not available if exact is set)")(
+      "weight-file,w", po::value<std::string>(&weight_file)->default_value(std::string()),
+      "Name of file containing a point weights. Format is one weight per line:\n  W1\n  ...\n  Wn ")(
       "output-file,o", po::value<std::string>(&output_file_diag)->default_value(std::string()),
       "Name of file in which the persistence diagram is written. Default print in std::clog")(
       "max-alpha-square-value,r", po::value<Filtration_value>(&alpha_square_max_value)
diff --git a/src/Alpha_complex/utilities/alphacomplex.md b/src/Alpha_complex/utilities/alphacomplex.md
index 527598a9..0d3c6027 100644
--- a/src/Alpha_complex/utilities/alphacomplex.md
+++ b/src/Alpha_complex/utilities/alphacomplex.md
@@ -46,6 +46,9 @@ for the Alpha complex construction.
 coefficient field Z/pZ for computing homology.

 * `-m [ --min-persistence ]` (default = 0) Minimal lifetime of homology feature

 to be recorded. Enter a negative value to see zero length intervals.

+* `-w [ --weight-file ]` is the path to the file containing the weights of the

+points (one value per line).

+Default version is not weighted.

 * `-e [ --exact ]` for the exact computation version.

 * `-f [ --fast ]` for the fast computation version.

 

@@ -58,6 +61,10 @@ to be recorded. Enter a negative value to see zero length intervals.
 N.B.:

 

 * Filtration values are alpha square values.

+* Weights values are explained on CGAL

+[dD Triangulations](https://doc.cgal.org/latest/Triangulation/index.html)

+and

+[Regular triangulation](https://doc.cgal.org/latest/Triangulation/index.html#title20) documentation.

 

 

 ## alpha_complex_3d_persistence ##