Added the files from the coxeter branch

4 files changed, 1001 insertions, 0 deletions

diff --git a/src/Coxeter_triangulation/include/gudhi/IO/Mesh_medit.h b/src/Coxeter_triangulation/include/gudhi/IO/Mesh_medit.h
new file mode 100644
index 00000000..b1cc0fe5
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/IO/Mesh_medit.h
@@ -0,0 +1,57 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Siargey Kachanovich
+ *
+ *    Copyright (C) 2019 Inria
+ *
+ *    Modification(s):
+ *      - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef IO_MESH_MEDIT_H_
+#define IO_MESH_MEDIT_H_
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+/* \class Mesh_medit
+ * \brief Structure to store a mesh that can be output in Medit .mesh file format
+ *  using the output_meshes_to_medit method.
+ *
+ * \ingroup coxeter_triangulation
+ */
+struct Mesh_medit {
+  /** \brief Type of a range of vertices. */
+  typedef std::vector<Eigen::VectorXd> Vertex_points;
+  /** \brief Type of a mesh element.
+   *  A pair consisting of a vector of vertex indices of type std::size_t
+   *  and of an integer that represents the common reference number for
+   *  the mesh elements of this type. */
+  typedef std::pair<std::vector<std::size_t>, std::size_t> Mesh_element;
+  /** \brief Type of a range of mesh elements. */
+  typedef std::vector<Mesh_element> Mesh_elements;
+  /** \brief Type of a range of scalar field . */
+  typedef std::vector<double> Scalar_field_range;
+
+  /** \brief Range of vertices of type Eigen::VectorXd to output. */
+  Vertex_points vertex_points;
+  /** \brief Range of edges. */
+  Mesh_elements edges;
+  /** \brief Range of triangles. */
+  Mesh_elements triangles;
+  /** \brief Range of tetrahedra. */
+  Mesh_elements tetrahedra;
+  /** \brief Range of scalar values over triangles. */
+  Scalar_field_range triangles_scalar_range;
+  /** \brief Range of scalar values over tetrahedra. */
+  Scalar_field_range tetrahedra_scalar_range;
+};
+
+} // namespace coxeter_triangulation 
+
+} // namespace Gudhi
+
+#endif
diff --git a/src/Coxeter_triangulation/include/gudhi/IO/build_mesh_from_cell_complex.h b/src/Coxeter_triangulation/include/gudhi/IO/build_mesh_from_cell_complex.h
new file mode 100644
index 00000000..5044660e
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/IO/build_mesh_from_cell_complex.h
@@ -0,0 +1,174 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Siargey Kachanovich
+ *
+ *    Copyright (C) 2019 Inria
+ *
+ *    Modification(s):
+ *      - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef IO_BUILD_MESH_FROM_CELL_COMPLEX_H_
+#define IO_BUILD_MESH_FROM_CELL_COMPLEX_H_
+
+#include <gudhi/IO/Mesh_medit.h>
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+struct Configuration {
+  bool toggle_edges = true,
+    toggle_triangles = true,
+    toggle_tetrahedra = true;
+  std::size_t ref_edges = 1,
+    ref_triangles = 1,
+    ref_tetrahedra = 1;
+
+  Configuration(bool t_edges, bool t_triangles, bool t_tetrahedra,
+		std::size_t r_edges, std::size_t r_triangles, std::size_t r_tetrahedra)
+    : toggle_edges(t_edges), toggle_triangles(t_triangles), toggle_tetrahedra(t_tetrahedra),
+      ref_edges(r_edges), ref_triangles(r_triangles), ref_tetrahedra(r_tetrahedra) {}
+};
+
+
+template <class Hasse_cell,
+	  class Simplex_cell_map>
+void populate_mesh(Mesh_medit& output,
+		   Simplex_cell_map& sc_map,
+		   Configuration configuration,
+		   std::size_t amb_d,
+		   std::map<Hasse_cell*, std::size_t> vi_map) {
+  using Mesh_element_vertices = Mesh_medit::Mesh_elements::value_type::first_type;
+  std::map<Hasse_cell*, std::size_t> ci_map;
+  std::size_t index = vi_map.size() + 1; // current size of output.vertex_points
+  if (sc_map.size() >= 3) 
+    for (const auto& sc_pair: sc_map[2]) {
+      Eigen::VectorXd barycenter = Eigen::VectorXd::Zero(amb_d);
+      std::set<std::size_t> vertex_indices;
+      Hasse_cell* cell = sc_pair.second;
+      for (const auto& ei_pair: cell->get_boundary())
+	for (const auto& vi_pair: ei_pair.first->get_boundary())
+	  vertex_indices.emplace(vi_map[vi_pair.first]);
+      for (const std::size_t& v: vertex_indices)
+	barycenter += output.vertex_points[v-1];
+      ci_map.emplace(std::make_pair(cell, index++));
+      output.vertex_points.emplace_back((1./vertex_indices.size()) * barycenter);    
+#ifdef GUDHI_COX_OUTPUT_TO_HTML
+      std::string vlist = " (" + std::to_string(index-1) + ")";
+      for (const std::size_t& v: vertex_indices)
+	vlist += " " + std::to_string(v);
+      cell_vlist_map.emplace(std::make_pair(to_string(cell), vlist));
+#endif
+    }
+
+  if (configuration.toggle_edges && sc_map.size() >= 2)
+    for (const auto& sc_map: sc_map[1]) {
+      Hasse_cell* edge_cell = sc_map.second;
+      Mesh_element_vertices edge;
+      for (const auto& vi_pair: edge_cell->get_boundary())
+	edge.push_back(vi_map[vi_pair.first]);
+      output.edges.emplace_back(std::make_pair(edge, configuration.ref_edges));
+#ifdef GUDHI_COX_OUTPUT_TO_HTML
+      std::string vlist;
+      for (const std::size_t& v: edge)
+	vlist += " " + std::to_string(v);
+      cell_vlist_map.emplace(std::make_pair(to_string(edge_cell), vlist));
+#endif
+    }
+  
+  if (configuration.toggle_triangles && sc_map.size() >= 3)
+    for (const auto& sc_pair: sc_map[2]) {
+      for (const auto& ei_pair: sc_pair.second->get_boundary()) {
+	Mesh_element_vertices triangle(1, ci_map[sc_pair.second]);
+	for (const auto& vi_pair: ei_pair.first->get_boundary())
+	  triangle.push_back(vi_map[vi_pair.first]);
+	output.triangles.emplace_back(std::make_pair(triangle, configuration.ref_triangles));
+      }
+    }
+  
+  if (configuration.toggle_tetrahedra && sc_map.size() >= 4)
+    for (const auto& sc_pair: sc_map[3]) {
+      Eigen::VectorXd barycenter = Eigen::VectorXd::Zero(amb_d);
+      std::set<std::size_t> vertex_indices;
+      Hasse_cell* cell = sc_pair.second;
+      for (const auto& ci_pair: cell->get_boundary())
+	for (const auto& ei_pair: ci_pair.first->get_boundary())
+	  for (const auto& vi_pair: ei_pair.first->get_boundary())
+	    vertex_indices.emplace(vi_map[vi_pair.first]);
+      for (const std::size_t& v: vertex_indices)
+	barycenter += output.vertex_points[v-1];
+      output.vertex_points.emplace_back((1./vertex_indices.size()) * barycenter);
+#ifdef GUDHI_COX_OUTPUT_TO_HTML
+      std::string vlist = " (" + std::to_string(index) + ")";
+      for (const std::size_t& v: vertex_indices)
+	vlist += " " + std::to_string(v);
+      cell_vlist_map.emplace(std::make_pair(to_string(cell), vlist));
+#endif
+
+      for (const auto& ci_pair: cell->get_boundary())
+	for (const auto& ei_pair: ci_pair.first->get_boundary()) {
+	  Mesh_element_vertices tetrahedron = {index, ci_map[sc_pair.second]};
+	  for (const auto& vi_pair: ei_pair.first->get_boundary())
+	    tetrahedron.push_back(vi_map[vi_pair.first]);
+	  output.tetrahedra.emplace_back(std::make_pair(tetrahedron, configuration.ref_tetrahedra));
+	}
+      index++;
+    }  
+}
+
+template <class Cell_complex>
+Mesh_medit build_mesh_from_cell_complex(const Cell_complex& cell_complex,
+					Configuration i_configuration = Configuration(),
+					Configuration b_configuration = Configuration()) {
+  using Hasse_cell = typename Cell_complex::Hasse_cell;
+  Mesh_medit output;
+  std::map<Hasse_cell*, std::size_t> vi_map; // one for vertices, other for 2d-cells
+  std::size_t index = 1; // current size of output.vertex_points
+
+  if (cell_complex.cell_point_map().empty())
+    return output;
+  std::size_t amb_d = std::min((int) cell_complex.cell_point_map().begin()->second.size(), 3);
+  
+  for (const auto& cp_pair: cell_complex.cell_point_map()) {
+#ifdef GUDHI_COX_OUTPUT_TO_HTML
+    std::string vlist;
+    vlist += " " + std::to_string(index);
+    cell_vlist_map.emplace(std::make_pair(to_string(cp_pair.first), vlist));
+#endif
+    vi_map.emplace(std::make_pair(cp_pair.first, index++));
+    output.vertex_points.push_back(cp_pair.second);
+    output.vertex_points.back().conservativeResize(amb_d);
+  }
+  
+
+  populate_mesh(output, cell_complex.interior_simplex_cell_maps(), i_configuration, amb_d, vi_map);
+#ifdef GUDHI_COX_OUTPUT_TO_HTML
+  for (const auto& sc_map: cell_complex.interior_simplex_cell_maps())    
+    for (const auto& sc_pair: sc_map) {
+      std::string simplex = "I" + to_string(sc_pair.first);
+      std::string cell = to_string(sc_pair.second);
+      std::string vlist = cell_vlist_map.at(cell).substr(1);
+      simplex_vlist_map.emplace(std::make_pair(simplex, vlist));
+    }
+#endif  
+  populate_mesh(output, cell_complex.boundary_simplex_cell_maps(), b_configuration, amb_d, vi_map);  
+#ifdef GUDHI_COX_OUTPUT_TO_HTML
+  for (const auto& sc_map: cell_complex.boundary_simplex_cell_maps())    
+    for (const auto& sc_pair: sc_map) {
+      std::string simplex = "B" + to_string(sc_pair.first);
+      std::string cell = to_string(sc_pair.second);
+      std::string vlist = cell_vlist_map.at(cell).substr(1);
+      simplex_vlist_map.emplace(std::make_pair(simplex, vlist));
+    }
+#endif  
+  return output;
+}
+
+} // namespace coxeter_triangulation 
+
+} // namespace Gudhi
+
+#endif
diff --git a/src/Coxeter_triangulation/include/gudhi/IO/output_debug_traces_to_html.h b/src/Coxeter_triangulation/include/gudhi/IO/output_debug_traces_to_html.h
new file mode 100644
index 00000000..43cca967
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/IO/output_debug_traces_to_html.h
@@ -0,0 +1,589 @@
+#ifdef GUDHI_DEBUG
+#define GUDHI_COX_OUTPUT_TO_HTML
+
+#include <sstream>
+#include <fstream>
+#include <vector>
+#include <list>
+#include <string>
+#include <regex>
+
+#include <Eigen/Dense>
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+template <class T>
+std::ostream& operator<<(std::ostream& os, const std::vector<T>& vector) {
+  os << "(";  
+  if (vector.empty()) {
+    os << ")";
+    return os;
+  }
+  auto v_it = vector.begin();
+  os << *v_it++;
+  for (; v_it != vector.end(); ++v_it)
+    os << ", " << *v_it;
+  os << ")";
+  return os;
+}
+
+/* A class to make the vector horizontal instead of vertical */
+struct Straighten {
+  Straighten (const Eigen::VectorXd& vector) : vector_(vector) {}
+  const Eigen::VectorXd& vector_;
+};
+
+std::ostream& operator<<(std::ostream& os, const Straighten& str) {
+  std::size_t size = str.vector_.size();
+  os << "(" << str.vector_(0);
+  if (size == 0) {
+    os << ")";
+    return os;
+  }
+  for (std::size_t i = 1; i < size; ++i)
+    os << ", " << str.vector_(i);
+  os << ")";
+  return os;
+}
+
+std::string id_from_simplex(const std::string& simplex) {
+  std::regex r("\\s+"), r2("\\(|\\)|\\{|\\}"), r3(","), r4("\\["), r5("\\]");
+  std::string output = std::regex_replace(simplex, r, "");
+  output = std::regex_replace(output, r2, ":");
+  output = std::regex_replace(output, r3, ".");
+  output = std::regex_replace(output, r4, "_");
+  output = std::regex_replace(output, r5, "");
+  return output;
+}
+
+template <typename T>
+std::string to_string(const T& t) {
+  std::ostringstream oss;
+  oss << t;
+  return oss.str();
+}
+
+struct MT_inserted_info {
+  std::string qr_face_, init_face_, qr_intersection_;
+  bool qr_success_, is_boundary_;
+  template <class Query_result,
+	    class Simplex_handle>
+  MT_inserted_info(const Query_result& qr, const Simplex_handle& face, bool is_boundary)
+    : qr_face_(to_string(face)), init_face_(to_string(face)),
+      qr_intersection_(to_string(qr.intersection)),
+      qr_success_(qr.success), is_boundary_(is_boundary) {}
+};
+std::list<MT_inserted_info> mt_seed_inserted_list, mt_inserted_list;
+
+struct CC_summary_info {
+  std::string face_, cell_;
+  template <class SC_pair>
+  CC_summary_info(const SC_pair& sc_pair)
+    : face_(to_string(sc_pair.first)), cell_(to_string(sc_pair.second)) {}
+};
+using CC_summary_list = std::list<CC_summary_info>;
+std::vector<CC_summary_list> cc_interior_summary_lists, cc_boundary_summary_lists;
+
+struct CC_detail_info {
+  enum class Result_type {self, face, coface, inserted, join_single, join_is_face};
+  std::string simplex_, trigger_, init_simplex_;
+  Result_type status_;
+  bool join_trigger_ = false;
+  std::list<std::string> faces_, post_faces_, cofaces_;
+  template <class Simplex_handle>
+  CC_detail_info(const Simplex_handle& simplex)
+    : simplex_(to_string(simplex)) {}
+};
+using CC_detail_list = std::list<CC_detail_info>;
+std::vector<CC_detail_list> cc_interior_detail_lists, cc_boundary_detail_lists;
+std::vector<CC_detail_list> cc_interior_insert_detail_lists, cc_boundary_insert_detail_lists;
+
+struct CC_prejoin_info {
+  enum class Result_type {join_single, join_is_face, join_different, join_same};
+  std::string simplex_, join_;
+  std::vector<std::string> faces_;
+  std::size_t dimension_;
+  Result_type status_;
+  template <class Simplex_handle>
+  CC_prejoin_info(const Simplex_handle& simplex)
+    : simplex_(to_string(simplex)), dimension_(simplex.dimension()) {}
+};
+using CC_prejoin_list = std::list<CC_prejoin_info>;
+std::vector<CC_prejoin_list> cc_interior_prejoin_lists, cc_boundary_prejoin_lists;
+
+
+struct CC_join_info {
+  enum class Result_type {self, face, coface, inserted, join_single, join_is_face};
+  std::string simplex_, join_, trigger_;
+  Result_type status_;
+  std::list<std::string> boundary_faces_;
+  std::list<std::string> faces_, post_faces_, cofaces_;
+  template <class Simplex_handle>
+  CC_join_info(const Simplex_handle& simplex)
+    : simplex_(to_string(simplex)) {}
+};
+bool join_switch = false;
+std::vector<CC_detail_list> cc_interior_join_detail_lists, cc_boundary_join_detail_lists;
+
+std::map<std::string, std::string> cell_vlist_map;
+std::map<std::string, std::string> simplex_vlist_map;
+
+std::ostringstream mt_ostream, vis_ostream;
+std::vector<std::ostringstream> cc_summary_ostream, cc_traces_ostream;
+
+std::string simplex_format(const std::string& simplex, bool is_boundary) {
+  std::string b_simplex = (is_boundary? "B": "I") + simplex;
+  std::string tooltiptext;
+  auto it = simplex_vlist_map.find(b_simplex);
+  if (it == simplex_vlist_map.end())
+    tooltiptext = "deleted";
+  else
+    tooltiptext = simplex_vlist_map.at(b_simplex);
+  return (std::string)"<a class=\"" + (is_boundary? "boundary": "interior")
+    + "\" href=\"#" + id_from_simplex(b_simplex) + "\">" + b_simplex
+    + "<span class=\"tooltiptext\">" + tooltiptext + "</span></a>";
+}
+
+std::string simplex_format(const std::string& b_simplex) {
+  bool is_boundary = b_simplex[0] == 'B';
+  std::string tooltiptext;
+  auto it = simplex_vlist_map.find(b_simplex);
+  if (it == simplex_vlist_map.end())
+    tooltiptext = "deleted";
+  else
+    tooltiptext = simplex_vlist_map.at(b_simplex);
+  return (std::string)"<a class=\"" + (is_boundary? "boundary": "interior")
+    + "\" href=\"#" + id_from_simplex(b_simplex) + "\">" + b_simplex
+    + "<span class=\"tooltiptext\">" + tooltiptext + "</span></a>";
+}
+
+
+void write_head(std::ofstream& ofs) {
+  ofs << "  <head>\n"
+      << "    <title>Cell complex debug trace</title>\n"
+      << "    <style>\n"
+      << "      a.boundary {\n"
+      << "        position: relative;\n"
+      << "        display: inline-block;\n"
+      << "	  color: darkred;\n"
+      << "	  background-color: lightgreen\n"
+      << "      }\n"
+      << "      a.interior {\n"
+      << "        position: relative;\n"
+      << "        display: inline-block;\n"
+      << "	  color: navy;\n"
+      << "	  background-color: yellow\n"
+      << "      }\n"
+      << "      .tooltiptext {\n"
+      << "        visibility: hidden;\n"
+      << "        width: 120px;\n"
+      << "        background-color: #555;\n"
+      << "        color: #fff;\n"
+      << "        text-align: center;\n"
+      << "        padding: 5px 0;\n"
+      << "        border-radius: 6px;\n"
+      << "        position: absolute;\n"
+      << "        z-index: 1;\n"
+      << "        bottom: 125%;\n"
+      << "        left: 50%;\n"
+      << "        margin-left: -60px;\n"
+      << "        opacity: 0;\n"
+      << "        transition: opacity 0.3s;\n"
+      << "      }\n"
+      << "      .boundary .tooltiptext::after {\n"
+      << "        content: \"\";\n"
+      << "        position: absolute;\n"
+      << "        top: 100%;\n"
+      << "        left: 50%;\n"
+      << "        margin-left: -5px;\n"
+      << "        border-width: 5px;\n"
+      << "        border-style: solid;\n"
+      << "        border-color: #555 transparent transparent transparent;\n"
+      << "      }\n"
+      << "      .interior .tooltiptext::after {\n"
+      << "        content: \"\";\n"
+      << "        position: absolute;\n"
+      << "        top: 100%;\n"
+      << "        left: 50%;\n"
+      << "        margin-left: -5px;\n"
+      << "        border-width: 5px;\n"
+      << "        border-style: solid;\n"
+      << "        border-color: #555 transparent transparent transparent;\n"
+      << "      }\n"
+      << "      .boundary:hover .tooltiptext {\n"
+      << "        visibility: visible;\n"
+      << "        opacity: 1;\n"
+      << "      }\n"    
+      << "      .interior:hover .tooltiptext {\n"
+      << "        visibility: visible;\n"
+      << "        opacity: 1;\n"
+      << "      }\n"    
+      << "      ul.nav {\n"
+      << "	  list-style-type: none;\n"
+      << "	  margin: 0;\n"
+      << "	  padding: 0;\n"
+      << "	  overflow: auto;\n"
+      << "	  background-color: #333;\n"
+      << "	  position: fixed;\n"
+      << "	  height: 100%;\n"
+      << "	  width: 15%;\n"
+      << "      }\n"
+      << "      ul.nav li a {\n"
+      << "	  display: block;\n"
+      << "	  color: white;\n"
+      << "	  text-align: left;\n"
+      << "	  padding: 14px 16px;\n"
+      << "	  text-decoration: none;\n"
+      << "      }\n"
+      << "      .active {\n"
+      << "	  background-color: #4CAF50;\n"
+      << "      }\n"
+      << "      div {\n"
+      << "        margin-left: 15%;\n"
+      << "        padding: 1px 16px\n"
+      << "      }\n"
+      << "      div.navi {\n"
+      << "        margin-left: 0%;\n"
+      << "        padding: 0px 0px\n"
+      << "      }\n"
+      << "      h1 {\n"
+      << "        margin-left: 15%;\n"
+      << "        padding: 1px 16px\n"
+      << "      }\n"
+      << "    </style>\n"
+      << "  </head>\n";
+}
+
+void write_nav(std::ofstream& ofs) {
+  ofs << "    <div class=\"navi\" style=\"margin-top:30px;background-color:#1abc9c;\">\n"
+      << "    <ul class=\"nav\">\n"
+      << "      <li><a href=\"#mant\">Manifold tracing</a></li>\n"
+      << "      <li><a href=\"#cell\">Cell complex</a>\n"
+      << "        <ul>\n";
+  for (std::size_t i = 0; i < cc_interior_summary_lists.size(); ++i) {
+    ofs << "          <li><a href=\"#dim" << i << "\">Dimension " << i << "</a>\n";
+    ofs << "            <ul>\n";
+    ofs << "              <li><a href=\"#dim" << i << "i\">Interior</a></li>\n";
+    if (i < cc_boundary_summary_lists.size()) {
+      ofs << "              <li><a href=\"#dim" << i << "b\">Boundary</a></li>\n";
+    }
+    ofs << "            </ul>\n";
+    ofs << "          </li>\n";
+  }
+  ofs << "        </ul>\n"
+      << "      </li>\n"
+      << "      <li><a href=\"#visu\">Visualization details</a></li>\n"
+      << "    </ul>\n"
+      << "    </div>\n";
+}
+
+void write_mt(std::ofstream& ofs) {
+  ofs << "    <div id=\"mant\">\n";
+  ofs << "      <h2> Manifold debug trace </h2>\n";
+  ofs << "      <h3> Simplices inserted during the seed phase </h3>\n";
+  ofs << "      <ul>\n";
+  for (const MT_inserted_info& mt_info: mt_seed_inserted_list) {
+    if (mt_info.qr_success_) {
+      ofs << "        <li>Inserted " << simplex_format(mt_info.qr_face_, mt_info.is_boundary_);
+      if (mt_info.qr_face_ != mt_info.init_face_)
+	ofs << " (initially " << simplex_format(mt_info.init_face_, mt_info.is_boundary_) << ")";
+      ofs << " intersection point is " << mt_info.qr_intersection_ << "</li>\n";
+    }
+    else
+      ofs << "        <li>Failed to insert " << mt_info.init_face_
+	  << "</li>\n";
+  }
+  ofs << "      </ul>\n";
+  ofs << "      <h3> Simplices inserted during the while loop phase </h3>\n";
+  ofs << "      <ul>\n";
+  for (const MT_inserted_info& mt_info: mt_inserted_list) {
+    if (mt_info.qr_success_) {
+      ofs << "        <li>Inserted " << simplex_format(mt_info.qr_face_, mt_info.is_boundary_);
+      if (mt_info.qr_face_ != mt_info.init_face_)
+	ofs << " (initially " << simplex_format(mt_info.init_face_, mt_info.is_boundary_) << ")";
+      ofs << " intersection point is " << mt_info.qr_intersection_ << "</li>\n";
+    }
+    else
+      ofs << "        <li>Failed to insert " << mt_info.init_face_
+	  << ")</li>\n";
+  }
+  ofs << "      </ul>\n";
+  ofs << "    </div>\n";
+}
+
+void write_cc(std::ofstream& ofs) {
+  ofs << "    <div id=\"cell\">\n"
+      << "      <h2> Cell complex debug trace </h2>\n"
+      << "      <p>Go to:</p>\n"
+      << "      <ul>\n";
+  for (std::size_t i = 0; i < cc_interior_summary_lists.size(); ++i) {
+    ofs << "	<li><a href=\"#dim" << i << "\">Dimension " << i << "</a></li>\n";
+  }
+  ofs << "      </ul>\n";
+  for (std::size_t i = 0; i < cc_interior_summary_lists.size(); ++i) {
+    ofs << "      <h3 id=\"dim" << i << "\"> Dimension " << i << "</h3>\n";
+    ofs << "      <h4 id=\"dim" << i << "i\"> Summary for interior simplices</h4>\n";
+    if (i < cc_boundary_summary_lists.size())
+      ofs << "      <p><a href=\"#dim" << i << "b\">Go to boundary</a></p>\n";
+    ofs << "        <ul>\n";
+    for (const CC_summary_info& cc_info: cc_interior_summary_lists[i])
+      ofs << "          <li id = \"" << id_from_simplex("I" + cc_info.face_) << "\">"
+	  << simplex_format(cc_info.face_, false)
+	  << " cell =" << cc_info.cell_ << "</li>\n";
+    ofs << "        </ul>\n";
+    ofs << "      <h4> Prejoin state of the interior cells of dimension " << i << "</h4>\n";
+    auto prejoin_it = cc_interior_prejoin_lists[i].begin();
+    while (prejoin_it != cc_interior_prejoin_lists[i].end()) {
+      std::size_t j = prejoin_it->dimension_;
+      ofs << "        <h5>" << j << "-dimensional ambient simplices</h5>\n";
+      ofs << "          <ul>\n";
+      for (; prejoin_it->dimension_ == j; ++prejoin_it) {
+	ofs << "          <li>" << simplex_format(prejoin_it->simplex_, false)
+	    << " join = " << simplex_format(prejoin_it->join_, false)
+	    << " boundary:\n"
+	    << "            <ul>\n";
+	for (const auto& face: prejoin_it->faces_)
+	  ofs << "              <li>" << simplex_format(face) << "</li>";
+	ofs << "            </ul>\n";
+	switch (prejoin_it->status_) {
+	case (CC_prejoin_info::Result_type::join_single):
+	  ofs << "            <p style=\"color: red\">Deleted "
+	      << simplex_format(prejoin_it->simplex_, false)
+	      << " as it has a single face.</p>";
+	  break;
+	case (CC_prejoin_info::Result_type::join_is_face):
+	  ofs << "            <p style=\"color: red\">Deleted "
+	      << simplex_format(prejoin_it->simplex_, false)
+	      << " as its join " << simplex_format(prejoin_it->join_, false)
+	      << " is one of the faces.</p>";
+	  break;
+	case (CC_prejoin_info::Result_type::join_different):
+	  ofs << "            <p style=\"color: magenta\">Deleted " << simplex_format(prejoin_it->simplex_, false)
+	      << " and replaced by its join " << simplex_format(prejoin_it->join_, false)
+	      << ".</p>";
+	  break;
+	case (CC_prejoin_info::Result_type::join_same):
+	  ofs << "            <p style=\"color: green\">Kept " << simplex_format(prejoin_it->simplex_, false)
+	      << ".</p>";
+	}
+	ofs << "          </li>";
+      }
+      ofs << "          </ul>\n";
+    }
+    ofs << "      <h4> Details for interior simplices</h4>\n";
+    ofs << "        <ul>\n";
+    for (const CC_detail_info& cc_info: cc_interior_detail_lists[i]) {
+      if (cc_info.status_ == CC_detail_info::Result_type::join_single) {
+	ofs << "          <li style=\"color:magenta\" id = \""
+	    << id_from_simplex("I" + cc_info.simplex_) << "\"> Simplex "
+	    << simplex_format(cc_info.simplex_, false) << " has only one face ("
+	    << simplex_format(cc_info.trigger_, false) << ") and is deleted.";
+	continue;
+      }
+      if (cc_info.status_ == CC_detail_info::Result_type::join_single) {
+	ofs << "          <li style=\"color:darkmagenta\" id = \""
+	    << id_from_simplex("I" + cc_info.simplex_) << "\"> The join of the simplex "
+	    << simplex_format(cc_info.simplex_, false) << " is one of its faces ("
+	    << simplex_format(cc_info.trigger_, false) << "), hence it is is deleted.";
+	continue;
+      }
+      ofs << "          <li> Insert_cell called for " << simplex_format(cc_info.simplex_, false)
+	  << "\n";
+      ofs << "            <ul>\n";
+      // for (const std::string& cof: cc_info.faces_)
+      // 	ofs << "              <li>Checking if " << simplex_format(cc_info.simplex_, false)
+      // 	    << " is a face of " << simplex_format(cof, false) << "\n";
+      ofs << "            </ul>\n";
+      ofs << "            <ul>\n";
+      if (cc_info.status_ == CC_detail_info::Result_type::self) {
+	ofs << "            <p><span style=\"color:blue\">The simplex "
+	    << simplex_format(cc_info.simplex_, false)
+	    << " already exists in the cell complex!</span></p>\n";
+      }
+      if (cc_info.status_ == CC_detail_info::Result_type::face) {
+	ofs << "            <p><span style=\"color:red\">The simplex "
+	    << simplex_format(cc_info.simplex_, false) << " is a face of the simplex "
+	    << simplex_format(cc_info.trigger_, false) << "!</span><br>\n";
+	ofs << "              <ul>\n";
+	for (const std::string post_face: cc_info.post_faces_)
+	  ofs << "                <li id = \"" << id_from_simplex("I" + post_face) << "\">"
+	      << "Post deleting " << simplex_format(post_face, false) << "</li>\n";
+	ofs << "              </ul>\n";
+	ofs << "            </p>\n";
+	ofs << "            <p id = \"" << id_from_simplex("I" + cc_info.trigger_) << "\">"
+	    << "Deleting " << simplex_format(cc_info.trigger_, false) << "</p>\n";
+      }
+      // for (const std::string& fac: cc_info.cofaces_)
+      // 	ofs << "              <li>Checking if " << simplex_format(cc_info.simplex_, false)
+      // 	    << " is a coface of " << simplex_format(fac, false) << "\n";
+      if (cc_info.status_ == CC_detail_info::Result_type::coface) {
+	ofs << "            <p><span style=\"color:darkorange\">The simplex "
+	    << simplex_format(cc_info.simplex_, false) << " is a coface of the simplex "
+	    << simplex_format(cc_info.trigger_, false) << "!</span><p>\n";
+      }
+      if (cc_info.status_ == CC_detail_info::Result_type::inserted) {
+	ofs << "            <p><span style=\"color:green\">Successfully inserted "
+	    << simplex_format(cc_info.simplex_, false) << "!</span><p>\n";
+      }
+      ofs << "            </ul>\n";
+      ofs << "          </li>\n";
+    }
+    ofs << "        </ul>\n";
+
+    if (i < cc_boundary_summary_lists.size()) {
+      ofs << "      <h4 id=\"dim" << i << "b\"> Summary for boundary simplices</h4>\n";
+      ofs << "      <p><a href=\"#dim" << i << "i\">Go to interior</a></p>\n";
+      ofs << "        <ul>\n";
+      for (const CC_summary_info& cc_info: cc_boundary_summary_lists[i])
+	ofs << "          <li  id = \"" << id_from_simplex("B" + cc_info.face_) << "\">"
+	    << simplex_format(cc_info.face_, true)
+	    << " cell =" << cc_info.cell_ << "</li>\n";
+      ofs << "        </ul>\n";
+      ofs << "      <h4> Prejoin state of the boundary cells of dimension " << i << "</h4>\n";
+      auto prejoin_it = cc_boundary_prejoin_lists[i].begin();
+      while (prejoin_it != cc_boundary_prejoin_lists[i].end()) {
+	std::size_t j = prejoin_it->dimension_;
+	ofs << "        <h5>" << j << "-dimensional ambient simplices</h5>\n";
+	ofs << "          <ul>\n";
+	for (; prejoin_it->dimension_ == j; ++prejoin_it) {
+	  ofs << "          <li>" << simplex_format(prejoin_it->simplex_, true)
+	      << " join = " << simplex_format(prejoin_it->join_, true)
+	      << " boundary:\n"
+	      << "            <ul>\n";
+	  for (const auto& face: prejoin_it->faces_)
+	    ofs << "              <li>" << simplex_format(face) << "</li>";
+	  ofs << "            </ul>\n";
+	  switch (prejoin_it->status_) {
+	  case (CC_prejoin_info::Result_type::join_single):
+	    ofs << "            <p style=\"color: red\">Deleted "
+		<< simplex_format(prejoin_it->simplex_, true)
+		<< " as it has a single face.</p>";
+	    break;
+	  case (CC_prejoin_info::Result_type::join_is_face):
+	    ofs << "            <p style=\"color: red\">Deleted "
+		<< simplex_format(prejoin_it->simplex_, true)
+		<< " as its join " << simplex_format(prejoin_it->join_, true)
+		<< " is one of the faces.</p>";
+	    break;
+	  case (CC_prejoin_info::Result_type::join_different):
+	    ofs << "            <p style=\"color: magenta\">Deleted " << simplex_format(prejoin_it->simplex_, true)
+		<< " and replaced by its join " << simplex_format(prejoin_it->join_, true)
+		<< ".</p>";
+	    break;
+	  case (CC_prejoin_info::Result_type::join_same):
+	    ofs << "            <p style=\"color: green\">Kept " << simplex_format(prejoin_it->simplex_, true)
+		<< ".</p>";
+	  }
+	  ofs << "          </li>";
+	}
+	ofs << "          </ul>\n";
+      }
+    }
+    if (i < cc_boundary_detail_lists.size()) {
+      ofs << "      <h4> Details for boundary simplices</h4>\n"
+	  << "        <ul>\n";
+      for (const CC_detail_info& cc_info: cc_boundary_detail_lists[i]) {
+	if (cc_info.status_ == CC_detail_info::Result_type::join_single) {
+	  ofs << "          <li style=\"color:magenta\" id = \""
+	      << id_from_simplex("B" + cc_info.simplex_) << "\"> Simplex "
+	      << simplex_format(cc_info.simplex_, true) << " has only one face ("
+	      << simplex_format(cc_info.trigger_, true) << ") and is deleted.";
+	  continue;
+	}
+	if (cc_info.status_ == CC_detail_info::Result_type::join_single) {
+	  ofs << "          <li style=\"color:darkmagenta\" id = \""
+	      << id_from_simplex("B" + cc_info.simplex_) << "\"> The join of the simplex "
+	      << simplex_format(cc_info.simplex_, true) << " is one of its faces ("
+	      << simplex_format(cc_info.trigger_, true) << "), hence it is is deleted.";
+	  continue;
+	}
+	ofs << "          <li> Insert_simplex called on " << simplex_format(cc_info.simplex_, true);
+	ofs << "            <ul>\n";
+	// for (const std::string& cof: cc_info.faces_)
+	//   ofs << "              <li>Checking if " << simplex_format(cc_info.simplex_, true)
+	//       << " is a face of " << simplex_format(cof, true) << "\n";
+	ofs << "            </ul>\n";
+	ofs << "            <ul>\n";
+	if (cc_info.status_ == CC_detail_info::Result_type::self) {
+	  ofs << "            <p><span style=\"color:blue\">The simplex "
+	      << simplex_format(cc_info.simplex_, true)
+	      << " already exists in the cell complex!</span></p>\n";
+	}
+	if (cc_info.status_ == CC_detail_info::Result_type::face) {
+	  ofs << "            <p><span style=\"color:red\">The simplex "
+	      << simplex_format(cc_info.simplex_, true) << " is a face of the simplex "
+	      << simplex_format(cc_info.trigger_, true) << "!</span><br>\n";
+	  ofs << "              <ul>\n";
+	  for (const std::string post_face: cc_info.post_faces_)
+	    ofs << "                <li id=\"" << id_from_simplex("B" + post_face)
+		<< "\">Post deleting " << simplex_format(post_face, true)
+		<< "</li>\n";
+	  ofs << "              </ul>\n";
+	  ofs << "            </p>\n";
+	  ofs << "            <p id=\"" << id_from_simplex(cc_info.trigger_)
+	      << "\">Deleting " << simplex_format(cc_info.trigger_, true) << "</p>\n";
+	}
+	// for (const std::string& fac: cc_info.cofaces_)
+	//   ofs << "              <li>Checking if " << simplex_format(cc_info.simplex_, true)
+	//       << " is a coface of " << simplex_format(fac, true) << "\n";
+	ofs << "            </ul>\n";
+	ofs << "          </li>\n";
+	if (cc_info.status_ == CC_detail_info::Result_type::coface) {
+	  ofs << "            <p><span style=\"color:darkorange\">The simplex "
+	      << simplex_format(cc_info.simplex_, true) << " is a coface of the simplex "
+	      << simplex_format(cc_info.trigger_, true) << "!</span><p>\n";
+	}
+	if (cc_info.status_ == CC_detail_info::Result_type::inserted) {
+	  ofs << "            <p><span style=\"color:green\">Successfully inserted "
+	      << simplex_format(cc_info.simplex_, true) << "!</span><p>\n";
+	}
+      }
+      ofs << "        </ul>\n";
+    }
+  }  
+  ofs << "    </div>\n";
+}
+
+void write_visu(std::ofstream& ofs) {
+  ofs << "    <div id=\"visu\">\n"
+      << "      <h2> Visualization details debug trace </h2>\n";
+  // std::vector<std::map<std::string, std::string> > vs_maps(cc_interior_summary_lists.size());
+  std::map<std::string, std::string> vs_map;
+  for (const auto& sv_pair: simplex_vlist_map)
+    vs_map.emplace(std::make_pair(sv_pair.second, sv_pair.first));
+  ofs << "      <ul>\n";
+  for (const auto& vs_pair: vs_map) {
+    std::string w_simplex = vs_pair.second.substr(1);
+    bool is_boundary = vs_pair.second[0] == 'B';
+    ofs << "      <li><b>" << vs_pair.first << "</b>: "
+	<< simplex_format(w_simplex, is_boundary) << "</li>\n";
+  }
+  ofs << "      </ul>\n";
+  ofs << "    </div>\n";
+}
+
+void write_to_html(std::string file_name) {
+  std::ofstream ofs(file_name + ".html", std::ofstream::out);
+  ofs << "<!DOCTYPE html>\n"
+      << "<html>\n";
+  write_head(ofs);
+  ofs << "  <body>\n";
+  write_nav(ofs);
+  ofs << "    <h1> Debug traces for " << file_name << " </h1>\n";
+  write_mt(ofs);
+  write_cc(ofs);
+  write_visu(ofs);  
+  ofs << "  </body>\n";
+  ofs << "</html>\n";
+
+  ofs.close();
+}
+
+} // namespace coxeter_triangulation
+
+} // namespace Gudhi
+
+#endif
+
diff --git a/src/Coxeter_triangulation/include/gudhi/IO/output_meshes_to_medit.h b/src/Coxeter_triangulation/include/gudhi/IO/output_meshes_to_medit.h
new file mode 100644
index 00000000..fe379242
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/IO/output_meshes_to_medit.h
@@ -0,0 +1,181 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Siargey Kachanovich
+ *
+ *    Copyright (C) 2019 Inria
+ *
+ *    Modification(s):
+ *      - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef IO_OUTPUT_MESHES_TO_MEDIT_H_
+#define IO_OUTPUT_MESHES_TO_MEDIT_H_
+
+#include <gudhi/IO/Mesh_medit.h>
+#include <Eigen/Dense>
+#include <fstream>
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+using Vertex_points = Mesh_medit::Vertex_points;
+using Mesh_elements = Mesh_medit::Mesh_elements;
+using Scalar_field_range = Mesh_medit::Scalar_field_range;
+
+template <std::size_t I = 0,
+	  typename... Meshes>
+typename std::enable_if<I == sizeof... (Meshes), void>::type
+fill_meshes (Vertex_points& vertex_points,
+	     Mesh_elements& edges,
+	     Mesh_elements& triangles,
+	     Mesh_elements& tetrahedra,
+	     Scalar_field_range& triangles_scalar_range,
+	     Scalar_field_range& tetrahedra_scalar_range,
+	     std::size_t index,
+	     const Meshes&... meshes) {
+}
+
+template <std::size_t I = 0,
+	  typename... Meshes>
+typename std::enable_if<I != sizeof... (Meshes), void>::type
+fill_meshes (Vertex_points& vertex_points,
+	     Mesh_elements& edges,
+	     Mesh_elements& triangles,
+	     Mesh_elements& tetrahedra,
+	     Scalar_field_range& triangles_scalar_range,
+	     Scalar_field_range& tetrahedra_scalar_range,
+	     std::size_t index,
+	     const Meshes&... meshes) {
+  auto mesh = std::get<I>(std::forward_as_tuple(meshes...));
+  for (const auto& v: mesh.vertex_points)
+    vertex_points.push_back(v);
+  for (const auto& e: mesh.edges) {
+    std::vector<std::size_t> edge;
+    for (const auto& v_i: e.first)
+      edge.push_back(v_i + index);
+    edges.emplace_back(std::make_pair(edge, e.second));
+  }
+  for (const auto& t: mesh.triangles) {
+    std::vector<std::size_t> triangle;
+    for (const auto& v_i: t.first)
+      triangle.push_back(v_i + index);
+    triangles.emplace_back(std::make_pair(triangle, t.second));
+  }
+  for (const auto& t: mesh.tetrahedra) {
+    std::vector<std::size_t> tetrahedron;
+    for (const auto& v_i: t.first)
+      tetrahedron.push_back(v_i + index);
+    tetrahedra.emplace_back(std::make_pair(tetrahedron, t.second));
+  }
+  for (const auto& b: mesh.triangles_scalar_range)
+    triangles_scalar_range.push_back(b);
+  for (const auto& b: mesh.tetrahedra_scalar_range)
+    tetrahedra_scalar_range.push_back(b);
+  fill_meshes<I+1, Meshes...>(vertex_points,
+			      edges,
+			      triangles,
+			      tetrahedra,
+			      triangles_scalar_range,
+			      tetrahedra_scalar_range,
+			      index + mesh.vertex_points.size(),
+			      meshes...);
+}
+
+/** \brief Outputs a text file with specified meshes that can be visualized in Medit.
+ *  
+ *  @param[in] amb_d Ambient dimension. Can be 2 or 3.
+ *  @param[in] file_name The name of the output file.
+ *  @param[in] meshes A pack of meshes to be specified separated by commas.
+ */
+template <typename... Meshes>
+void output_meshes_to_medit(std::size_t amb_d, std::string file_name, const Meshes&... meshes) {
+  Vertex_points vertex_points;
+  Mesh_elements edges, triangles, tetrahedra;
+  Scalar_field_range triangles_scalar_range, tetrahedra_scalar_range;
+  fill_meshes(vertex_points,
+	      edges,
+	      triangles,
+	      tetrahedra,
+	      triangles_scalar_range,
+	      tetrahedra_scalar_range,
+	      0,
+	      meshes...);
+  
+  std::ofstream ofs (file_name + ".mesh", std::ofstream::out);
+  std::ofstream ofs_bb (file_name + ".bb", std::ofstream::out);
+  
+  if (amb_d == 2) {
+    ofs << "MeshVersionFormatted 1\nDimension 2\n";
+    ofs_bb << "2 1 ";
+    ofs << "Vertices\n" << vertex_points.size() << "\n";
+    for (auto p: vertex_points) {
+      ofs << p[0] << " " << p[1] << " 2\n";
+    }
+    ofs << "Edges " << edges.size() << "\n";
+    for (auto e: edges) {
+      for (auto v: e.first)      
+    	ofs << v << " ";
+      ofs << e.second << std::endl;
+    }
+    ofs << "Triangles " << triangles.size() << "\n";
+    for (auto s: triangles) {
+      for (auto v: s.first) {
+    	ofs << v << " ";
+      }
+      ofs << s.second << std::endl;
+    }
+
+    ofs_bb << triangles_scalar_range.size() << " 1\n";
+    for (auto& b: triangles_scalar_range)
+      ofs_bb << b << "\n";
+
+  }
+  else {
+    ofs << "MeshVersionFormatted 1\nDimension 3\n";
+    ofs_bb << "3 1 ";
+    ofs << "Vertices\n" << vertex_points.size() << "\n";
+    for (auto p: vertex_points) {
+      ofs << p[0] << " " << p[1] << " " << p[2] << " 2\n";
+    }
+    ofs << "Edges " << edges.size() << "\n";
+    for (auto e: edges) {
+      for (auto v: e.first)      
+    	ofs << v << " ";
+      ofs << e.second << std::endl;
+    }
+    ofs << "Triangles " << triangles.size() << "\n";
+    for (auto s: triangles) {
+      for (auto v: s.first) {
+    	ofs << v << " ";
+      }
+      ofs << s.second << std::endl;
+    }
+    ofs << "Tetrahedra " << tetrahedra.size() << "\n";
+    for (auto s: tetrahedra) {
+      for (auto v: s.first) {
+    	ofs << v << " ";
+      }
+      ofs << s.second << std::endl;
+    }
+
+    ofs_bb << triangles_scalar_range.size() + tetrahedra_scalar_range.size() << " 1\n";
+    for (auto& b: triangles_scalar_range)
+      ofs_bb << b << "\n";
+    for (auto& b: tetrahedra_scalar_range)
+      ofs_bb << b << "\n";
+  }
+
+  
+  ofs.close();
+  ofs_bb.close();
+
+}
+
+} // namespace coxeter_triangulation 
+
+} // namespace Gudhi
+
+#endif