Merge 179f2dfaf2 into 26a5fc70e4

.gitignore

@@ -1172,3 +1172,5 @@ patches_after_merge.ply
 CMakeUserPresets.json
 /.cache
 compile_commands.json
+/HDVF/test/HDVF/tmp
+/HDVF/examples/HDVF/tmp

Algebraic_foundations/doc/Algebraic_foundations/CGAL/number_utils.h

@@ -22,7 +22,7 @@ namespace CGAL {
 \ingroup PkgAlgebraicFoundationsRef
 
 The template function `compare()` compares the first argument with respect to
-the second, i.e.\ it returns `CGAL::LARGER` if \f$ x\f$ is larger than \f$ y\f$.
+the second.
 
 In case the argument types `NT1` and `NT2` differ,
 `compare` is performed with the semantic of the type determined via

Documentation/doc/Documentation/packages.txt

@@ -131,7 +131,11 @@
 \package_listing{Stream_lines_2}
 \package_listing{Classification}
 \package_listing{Heat_method_3}
+
+\cgalPackageSection{PartTopology,Topology}
+
 \package_listing{Surface_mesh_topology}
+\package_listing{HDVF}
 
 \cgalPackageSection{PartSearchStructures,Spatial Searching and Sorting}
 

HDVF/doc/HDVF/Concepts/AbstractChainComplex.h

@@ -0,0 +1,136 @@
+/*!
+\ingroup PkgHDVFConcepts
+\cgalConcept
+
+The concept `AbstractChainComplex` describes the requirements for (topological) chain complexes associated to abstract complexes used in the concept `HDVF`.
+
+It provides methods to:
+
+ - get the dimension of the complex, the number of cells in each dimension
+ - get the boundary and co-boundary of cell(s)
+ - get the vertices of a given cell
+ - output the complex in text format
+
+ Cells are indexed along each dimension and thus identified by their index together with their dimension.
+
+
+\cgalHasModelsBegin
+\cgalHasModelsBare{`CGAL::Homological_discrete_vector_field::Abstract_simplicial_chain_complex<IntegralDomainWithoutDivision>`}
+\cgalHasModelsBare{`CGAL::Homological_discrete_vector_field::Sub_chain_complex_mask<IntegralDomainWithoutDivision, AbstractChainComplex>`}
+\cgalHasModelsEnd
+
+*/
+
+class AbstractChainComplex
+{
+public:
+/// \name Types
+/// @{
+
+/*!
+Type of coefficients ring used to compute homology, model of `IntegralDomainWithoutDivision`
+ */
+typedef CoefficientRing Coefficient_ring;
+
+/*!
+ Type of column-major chains (returned by the boundary operator)
+ */
+typedef CGAL::OSM::Sparse_chain<CoefficientRing, CGAL::OSM::COLUMN> Column_chain;
+
+/*!
+ Type of row-major chains (returned by the co-boundary operator)
+ */
+typedef CGAL::OSM::Sparse_chain<CoefficientRing, CGAL::OSM::ROW> Row_chain ;
+
+/*!
+ Type of column-major sparse matrices (used to store the boundary operator)
+ */
+typedef CGAL::OSM::Sparse_matrix<CoefficientRing, CGAL::OSM::COLUMN> Column_matrix;
+
+/// @}
+
+/// \name Operators
+/// @{
+
+/*!
+Assignment operator.
+
+The operator creates a copy of `complex`.
+ */
+    AbstractChainComplex& operator= (const AbstractChainComplex& complex);
+
+/// @}
+
+/// \name Access functions
+/// @{
+
+/*!
+Returns the dimension of the complex, that is, the largest dimension of cells.
+ */
+int dimension();
+
+/*!
+Returns the number of cells of dimension `q`.
+If `q` is negative of larger than the dimension of the complex, returns 0.
+ */
+size_t number_of_cells(int q);
+
+/*!
+\brief Returns all boundary matrices.
+
+The function returns constant reference to a vector of column-major sparse matrices. The `q`-th element of this vector is the matrix of \f$\partial_q\f$, which gives the boundary of cells of dimension `q`(as a linear combination of `q`-1 cells).
+
+ */
+const vector<Column_matrix> & boundary_matrices() const;
+
+/*!
+ \brief Returns the boundary matrix of dimension `q` (ie.\ the matrix of \f$\partial_q\f$).
+
+The function returns a column-major sparse matrices.
+ */
+const Column_matrix & boundary_matrix(int q) const;
+
+/*!
+ \brief Returns the boundary of the cell of index `id_cell` in dimension `q`.
+
+This boundary is a finite linear combination of cells of dimension `q`-1. It is encoded as a column-major chain (which maps each cell with a non-zero coefficient to this coefficient). This boundary  is thus the `id_cell`-th column of the boundary matrix in dimension `q`.
+
+*/
+Column_chain d(size_t id_cell, int q);
+
+
+/*!
+ \brief Returns the co-boundary of the cell of index `id_cell` in dimension `q`.
+
+This boundary is a finite linear combination of cells of dimension `q`+1. It is encoded as a row-major chain (which maps each cell with a non-zero coefficient to this coefficient). This co-boundary  is thus the `id_cell`-th row of the boundary matrix in dimension `q`+1.
+*/
+Row_chain cod(size_t id_cell, int q);
+
+/*!
+ \brief Returns the vertices of a given cell (that is, the indices of its faces of dimension 0).
+
+*/
+std::vector<size_t> bottom_faces(size_t id_cell, int q) const;
+
+/*!
+ * \brief Returns the cofaces of a given chain in dimension `q`.
+ *
+ * The resulting chain, whatever the storage format of the input, is column-major, lies in dimension `q`+1 and is null if this dimension exceeds the dimension of the complex.
+ *
+ * \tparam CoefficientR `CoefficientRing` of the chain.
+ * \tparam StorageF `StorageFormat` of the chain.
+*/
+template <typename CoefficientR, int StorageF>
+Column_chain cofaces_chain (SparseChain<CoefficientR, StorageF> chain, int q) const;
+
+/// @}
+
+};
+
+
+/*!
+\relates AbstractChainComplex
+Inserts the chain complex in text mode in the stream.
+*/
+std::ostream& operator<<(std::ostream& out, const AbstractChainComplex& complex) const;
+

HDVF/doc/HDVF/Concepts/Filtration.h

@@ -0,0 +1,181 @@
+/*!
+\ingroup PkgHDVFConcepts
+\cgalConcept
+
+The concept `Filtration` describes the requirements for persistent filtrations associated to persistent homology computation.
+
+A filtration is associated to an `AbstractChainComplex`. Each cell is equiped with a scalar value (called its degree) and the filtration is an enumeration of cells in any dimension by increasing degrees.
+
+A filtration class provides:
+
+ - an iterator to visit all cells by increasing degrees.
+ - getters to get the degree of a cell, and the cell of a given index along the filtration.
+ - an overload of `operator<<()` to output filtrations.
+
+Cells are indexed along each dimension and thus identified by their index together with their dimension.
+
+\cgalHasModelsBegin
+\cgalHasModelsBare{`CGAL::Homological_discrete_vector_field::Filtration_core<ChainComplex,Degree>`}
+\cgalHasModelsBare{`CGAL::Homological_discrete_vector_field::Filtration_lower_star<ChainComplex,Degree>`}
+\cgalHasModelsEnd
+
+*/
+
+class Filtration
+{
+public:
+    /*! \brief Type of the underlying complex. */
+    typedef unspecified_type Complex;
+
+    /*! \brief (Scalar) type of degrees (model of `RealEmbeddable`).
+     */
+    typedef unspecified_type Degree ;
+
+    /*! \brief Type of coefficients used to compute homology. */
+    typedef Complex::Coefficient_ring Coefficient_ring;
+
+    /*! \brief Type for indexing uniquely a cell.
+     *
+     * As stated in `AbstractChainComplex`, cells are identified by their dimension, together with their index along this dimension. The type `Cell` stores this pair:
+     *
+     * - First element of the pair: index of the cell.
+     * - Second element of the pair: dimension of the cell.
+     */
+    typedef std::pair<size_t, int> Cell ;
+
+    /*! \brief Value type of the filtration iterator.
+     * Contains a cell (identified by its index and dimension in a `Cell`) and its associated degree.
+     */
+    typedef struct {
+        Cell cell_dim ;
+        Degree degree ;
+    } Filtration_iter_value ;
+
+protected:
+    /*!
+     Type of column-major sparse matrices.
+     */
+    typedef SparseMatrix<Coefficient_ring,CGAL::OSM::COLUMN> Column_matrix ;
+
+    /*!
+     Type of row-major sparse matrices.
+     */
+    typedef SparseMatrix<Coefficient_ring,CGAL::OSM::ROW> Row_matrix ;
+
+    /*!
+     Type of column-major chains.
+     */
+    typedef SparseChain<Coefficient_ring,CGAL::OSM::COLUMN> Column_chain ;
+
+    /*!
+     Type of row-major chains.
+     */
+    typedef SparseChain<Coefficient_ring,CGAL::OSM::ROW> Row_chain ;
+
+public:
+    /**
+     * \brief The iterator over filtration.
+     *
+     * Iterate the filtration by increasing degrees.
+     */
+    struct iterator
+    {
+        // Iterator tags
+        using iterator_category = std::forward_iterator_tag;
+        using difference_type   = std::ptrdiff_t;
+        using value_type        = Filtration_iter_value;
+
+        /*! Iterator constructor
+         */
+        iterator(const Filtration& f, size_t i=0);
+
+        /*! Dereference operator.
+         */
+        value_type operator*() const;
+
+        /*! Gets the index (time) associated to the iterator.
+         */
+        size_t time () const ;
+
+        /*! Gets the `Cell` (cell and its dimension) associated to the iterator.
+         */
+        Cell cell_dim () const ;
+
+        /*! Gets the degree associated to the iterator.
+         */
+        Degree degree () const ;
+
+        /*!
+         * \brief Prefix incrementation.
+         * Move to next index in the filtration.
+         */
+        iterator& operator++();
+
+        /*! \brief Postfix incrementation.
+         */
+        iterator operator++(int);
+
+        /*!
+         * \brief Equality check.
+         * \returns True if the indices are equal.
+         */
+        friend bool operator== (const iterator& a, const iterator& b);
+
+        /*!
+         * \brief Inequality check.
+         * \returns True if the indices are different.
+         */
+        friend bool operator!= (const iterator& a, const iterator& b);
+    };
+
+    /*!
+     * \brief Iterator to the beginning of the filtration.
+     */
+    iterator begin();
+
+    /*!
+     * \brief Returns a past-the-end iterator.
+     */
+    iterator end();
+
+    // getters
+    /*! \brief Gets the filtration size.
+     */
+    size_t size () const;
+
+    /*! \brief Gets the cell (that is cell index and dimension) at the index `i` of the filtration.
+     */
+    Cell cell_index_dimension (size_t i) const;
+
+    /*! \brief Gets the degree of the `i`th element of the filtration.
+     */
+    Degree degree (size_t i) const;
+
+    // Filtration verification
+    /*! \brief Checks that a filtration is valid.
+     * Checks that cells are ordered in increasing degrees and all cells have indices larger than their faces.
+     */
+    bool is_valid() const;
+
+    // Input/output filtration
+
+
+    /**
+     * \brief Exports the filtration time indices.
+     *
+     * The method exports the time index of every cell in each dimension.
+     */
+    vector<vector<size_t> > export_filtration () const
+};
+
+    /*!
+    \relates Filtration
+    \brief extracts a filtration from a stream.
+     */
+    istream & operator>>(istream & in, Filtration &f);
+
+    /*!
+    \relates Filtration
+    \brief Inserts a filtration into a stream.
+     */
+    ostream & operator<<(ostream & out, const Filtration &f);

HDVF/doc/HDVF/Concepts/GeometricChainComplex.h

@@ -0,0 +1,45 @@
+/*!
+\ingroup PkgHDVFConcepts
+\cgalConcept
+
+The concept `GeometricChainComplex` refines the concept `AbstractChainComplex` and describes the requirements for (topological) chain complexes associated to geometric complexes used in the concept `HDVF`.
+ It adds to `AbstractChainComplex` methods to get vertex coordinates.
+
+ \cgalRefines{AbstractChainComplex}
+
+\cgalHasModelsBegin
+\cgalHasModelsBare{`CGAL::Homological_discrete_vector_field::Simplicial_chain_complex<CoefficientRing,Traits>`}
+\cgalHasModelsBare{`CGAL::Homological_discrete_vector_field::Cubical_chain_complex<CoefficientRing,Traits>`}
+\cgalHasModelsEnd
+
+*/
+
+class GeometricChainComplex
+{
+public:
+/// \name Types
+/// @{
+
+/*!
+ Type of coordinates (the vector size is free, hence coordinates can be any dimension).
+ */
+typedef unspecified_type Point ;
+/// @}
+
+
+/// \name Access functions
+/// @{
+
+/*!
+Returns the vector of vertex coordinates.
+ */
+const std::vector<Point>& points() const;
+
+/*!
+ Returns the coordinates of the cell of index `i` and dimension 0.
+ */
+Point point(size_t i) const;
+
+/// @}
+
+};

HDVF/doc/HDVF/Concepts/HDVF.h

@@ -0,0 +1,414 @@
+/*!
+\ingroup PkgHDVFConcepts
+\cgalConcept
+
+The concept `HDVF` describes the requirements for Homological Discrete Vector Fields (%HDVF for short) , a theory of computational homology unifying discrete Morse theory and effective homology. HDVFs were introduced by Aldo Gonzalez-Lorenzo in his PhD (see [AGL,2017], [AGL,2016]).
+
+A HDVF is a combinatorial tool computing the homology of a chain complex (described by the `AbstractChainComplex` concept) including "geometric" information such as homology and cohomology generators.
+
+ As described in the user manual description, an HDVF is:
+ - a partition of the cells of the underlying complex into three classes: primary (P), secondary (S) and critical (C), together with partial invertibility conditions for the boundary operator,
+ - an associated reduction \f$(f,g,h)\f$ from the chain complex onto a chain complex generated by critical cells and a reduced boundary.
+ A HDVF is perfect when the reduced boundary is null; in this case, the maps \f$f\f$ and \f$g\f$ (and their duals) provide homology and cohomology information (generators, annotation...).
+
+ More precisely, a HDVF over a chain complex C, derived from a complex K, is a partition of the cells: \f$K = P\coprod S\coprod C\f$ with \f$P\f$ primary cells, \f$S\f$ secondary cells and \f$C\f$ critical cells, such that the sub-matrix:
+ \f\[\partial (S)|_P\text{ is invertible.}
+ \f\]
+ In what follows, given \f$q\in \mathbb N\f$ and \f$G\f$ a sub-group of the qth chain group \f$C_q\f$, we denote by:
+
+ - \f$i_G\f$ the inclusion \f$G\to C_q\f$
+ - \f$j_G\f$ the projection \f$C_q\to G\f$
+
+That is, for all \f$q\in \mathbb N\f$, if
+
+ \f\[\mathrm{Mat}(\partial_q) = \begin{array}{lc}
+ & P_q \ \ \ \ S_q \ \ \ \ C_q \\
+ \begin{array}{l}
+ P_{q-1} \\
+ S_{q-1} \\
+ C_{q-1} \\
+ \end{array} & \left(\begin{array}{ccc}
+ \star & D_{PS} & \star \\
+\star & \star & \star \\
+\star & \star & \star \\
+ \end{array}\right)\\
+ \end{array}
+ \f\]
+ then \f$D_{PS}\f$, the matrix of \f$j_P\circ \partial\circ i_S\f$, also written \f$\partial (S)|_P\f$ is invertible.
+
+ In order to build a perfect HDVF, a pairing operation called `A()` associates (under conditions) two critical cells (becoming PRIMARY and SECONDARY) and updates the reduction in time \f$\mathscr O(n^2)\f$. Intuitively, this operation bears some similitudes (see users guide) with discrete Morse theory's DGVF arrows (and an HDVF can always be represented by such a vector field - possibly with cycles).
+ Given a HDVF \f$X = (P,S,C)\f$, A operation between two cells \f$(\gamma_1, \gamma_2)\f$ of respective dimension \f$q/q+1\f$ is valid if \f$(P',S',C') = (P\cup\{\gamma_1\}, P\cup\{\gamma_2\}, C\backslash\{\gamma_1,\gamma_2\})\f$ is a HDVF (ie.\ if \f$j_{P'}\circ \partial\circ i_{S'}\f$ is still invertible). This is actually the case if and only if \f$\langle d(\gamma_2), \gamma_1 \rangle\f$ (with \f$d\f$ the reduced boundary operator) is invertible. This condition is called the *validity condition for A*; all the `find_pair_A()` and `find_pairs_A()` methods search for such valid pairs.
+
+
+ Starting from this operation, `HDVF` provides two methods to generate perfect HDVFs (and hence compute homology):
+ - `compute_perfect_hdvf()` which computes a perfect HDVF by chosing at each step the first A-pairings available (depending on cells ordering in the chain complex)
+ - `compute_rand_perfect_hdvf()` which computes a perfect random HDVF by chosing at each step a rand A-pairing among all possible ones  (this option is thus slower)
+
+ For efficiency, the HDVF reduction can be computed partially:
+ - compute only the reduced boundary  (option `OPT_BND`): to get Betti numbers
+ - compute \f$\partial'\f$ together with \f$h\f$ and \f$g\f$ (option `OPT_G`): to get homology generators
+ - compute \f$\partial'\f$ together with \f$h\f$ and \f$f\f$ (option `OPT_F`): to get co-homology generators
+ - compute the full reduction (option `OPT_FULL`): to get the full homological / cohomological information
+
+ Let us consider a simple cubical complex example (see \ref secHDVFcubical_chain_complex for an introduction). In what follows, PRIMARY cells are coloured green, SECONDARY are coloured red and critical blue. The picture below shows an initial empty HDVF where all cells are critical:
+
+ <img src="HDVF_ex1.png" align="center" width=40%/>
+
+ The boundary matrix of \f$\partial_2\f$ (which equals the reduced boundary matrix as the HDVF is empty) for \f$\mathbb Z\f$-homology is (cells are identified by their Khalimsky coordinates):
+ \f\[\mathrm{Mat}(\partial_2) = \begin{array}{lc}
+ & (3,1)\\
+ \begin{array}{l}
+ (3,0) \\
+ (2,1) \\
+ (4,1) \\
+ (3,2) \\
+ \end{array} & \left(\begin{array}{c}
+ 1 \\
+ -1 \\
+ 1 \\
+ -1 \\
+ \end{array}\right)\\
+ \end{array}
+ \f\]
+ Let us consider \f$\gamma_2 = (3,1)\f$ and \f$\gamma_1 = (3,2)\f$, the condition for A is valid (-1 is invertible). The following images illustrates this pairing in discrete Morse theory "style" (pink arrow in left complex) and the resulting HDVF (right complex).
+
+ <img src="HDVF_ex1_A1_arrow.png" align="center" width=30%/>
+ <img src="HDVF_ex1_A1.png" align="center" width=30%/>
+
+ The reduced boundary matrix in dimension 2 becomes empty and the matrix of \f$d_1\f$ is similar to the boundary matrix \f$\partial_1\f$ with the row/colums associated to previous cells removed.
+
+ Thus, if we now consider \f$\gamma_2 = (4,1)\f$ and \f$\gamma_1 = (4,0)\f$, the condition for A is valid; cells can be paired and we obtain the following HDVF:
+
+ <img src="HDVF_ex1_A2_arrow.png" align="center" width=30%/>
+ <img src="HDVF_ex1_A2.png" align="center" width=30%/>
+
+ The process can be continued step by step:
+
+ <img src="HDVF_ex1_A3_arrow.png" align="center" width=30%/>
+ <img src="HDVF_ex1_A3.png" align="center" width=30%/>
+
+ Until we eventually get a perfect HDVF:
+
+ <img src="HDVF_ex1_perfect.png" align="center" width=40%/>
+
+ Each hole is identified by a critical cell and the \f$g\f$ map of the reduction provides associated homology generators (highlighted in blue):
+
+ <img src="HDVF_ex1_perfect_hom1.png" align="center" width=30%/>
+ <img src="HDVF_ex1_perfect_hom2.png" align="center" width=30%/>
+ <img src="HDVF_ex1_perfect_hom3.png" align="center" width=30%/>
+
+ A corresponding cohomology generators (highlighted in pink):
+
+ <img src="HDVF_ex1_perfect_cohom1.png" align="center" width=30%/>
+ <img src="HDVF_ex1_perfect_cohom2.png" align="center" width=30%/>
+ <img src="HDVF_ex1_perfect_cohom3.png" align="center" width=30%/>
+
+
+\cgalHasModelsBegin
+\cgalHasModelsBare{`CGAL::Homological_discrete_vector_field::Hdvf_core<IntegralDomainWithoutDivision, AbstractChainComplex, SparseChain, Sparse_matrix>`}
+\cgalHasModelsBare{`CGAL::Homological_discrete_vector_field::Hdvf<IntegralDomainWithoutDivision, AbstractChainComplex>`}
+\cgalHasModelsBare{`CGAL::Homological_discrete_vector_field::Hdvf_persistence<IntegralDomainWithoutDivision, AbstractChainComplex>`}
+\cgalHasModelsBare{`CGAL::Homological_discrete_vector_field::Hdvf_duality<IntegralDomainWithoutDivision, AbstractChainComplex>`}
+\cgalHasModelsEnd
+
+ \sa `AbstractChainComplex`
+ \sa `IntegralDomainWithoutDivision`
+ \sa `SparseChain`
+ \sa `SparseMatrix`
+
+ *How to describe constants declared in the namespace Homological_discrete_vector_field and used everywhere? PSC_flag, options, exporttype*
+
+ [AGL, 2017] Aldo Gonzalez-Lorenzo, Alexandra Bac, Jean-Luc Mari, Pedro Real. Allowing cycles in discrete Morse theory, Topology and its Applications, Volume 228, 2017, Pages 1-35.
+
+ [AGL, 2016] Aldo Gonzalez-Lorenzo. Computational Homology Applied to Discrete Objects. Discrete Mathematics [cs.DM]. Aix-Marseille Université; Universidad de Sevilla, 2016.
+
+*/
+
+class HDVF
+{
+public:
+/// \name Types
+/// @{
+
+/*!
+ * \brief Type of underlying chain complex (a model of `AbstractChainComplex`).
+ */
+typedef unspecified_type Complex;
+
+/*! \brief Type of coefficients used to compute homology. */
+typedef Complex::Coefficient_ring Coefficient_ring;
+
+/*!
+ * \brief Type of sparse chains (a model of `SparseChain`).
+ */
+    typedef unspecified_type Sparse_chain;
+
+/*!
+ * \brief Type of sparse matrices (a model of `SparseMatrix`).
+ */
+    typedef unspecified_type Sparse_matrix;
+
+/*! \brief Structure to represent data for HDVF operations (pairs of cells).
+ *
+ * Cells are always sorted so that the dimension of `sigma` is lesser than the dimension of `tau`.
+ */
+struct Cell_pair {
+    /// Index of the first cell
+    size_t sigma;
+    /// Index of the second cell
+    size_t tau;
+    /// Dimension of cells: `dim`/`dim`+1 for `A()` and `R()`, `dim`/`dim` for other operations
+    int dim;
+};
+
+/*!
+ Type of column-major chains
+ */
+typedef Sparse_chain<CoefficientType, CGAL::OSM::COLUMN> Column_chain;
+
+/*!
+ Type of row-major chains
+ */
+typedef Sparse_chain<CoefficientType, CGAL::OSM::ROW> Row_chain ;
+
+/*!
+ Type of column-major sparse matrices
+ */
+typedef Sparse_matrix<CoefficientType, CGAL::OSM::COLUMN> Column_matrix;
+
+/*!
+ Type of row-major sparse matrices
+ */
+typedef Sparse_matrix<CoefficientType, CGAL::OSM::ROW> Row_matrix;
+
+/// @}
+
+/// \name Functions to find valid pairs for `A()` operations
+/// @{
+
+/*!
+ * \brief Finds a valid Cell_pair of dimension q / q+1 for A.
+ *
+ * The function searches a pair of critical cells \f$(\gamma_1, \gamma2)\f$ of dimension q / q+1, valid for A (ie.\ such that \f$\langle \partial_{q+1}(\gamma_2), \gamma_1 \rangle\f$ invertible). It returns the first valid pair found by iterators.
+ *
+ * \param q Lower dimension of the pair.
+ * \param found Reference to a %Boolean variable. The method sets `found` to `true` if a valid pair is found, `false` otherwise.
+ */
+virtual Cell_pair find_pair_A(int q, bool &found) const;
+
+/*!
+ * \brief Finds a valid Cell_pair for A containing `gamma` (a cell of dimension `q`)
+ *
+ * The function searches a cell \f$\gamma'\f$ such that one of the following conditions holds:
+ * - \f$\gamma'\f$ has dimension q+1 and \f$(\gamma, \gamma')\f$ is valid for A (ie.\ such that \f$\langle \partial_{q+1}(\gamma'), \gamma \rangle\f$ invertible),
+ * - \f$\gamma'\f$ has dimension q-1 and \f$(\gamma', \gamma)\f$ is valid for A (ie.\ such that \f$\langle \partial_{q}(\gamma), \gamma' \rangle\f$ invertible).
+ *
+ * \param q Dimension of the cell `gamma`.
+ * \param found Reference to a %Boolean variable. The method sets `found` to `true` if a valid pair is found, `false` otherwise.
+ * \param gamma Index of a cell to pair.
+ */
+virtual Cell_pair find_pair_A(int q, bool &found, size_t gamma) const;
+
+/**
+ * \brief Finds *all* valid Cell_pair of dimension q / q+1 for A.
+ *
+ * The function searches all pairs of critical cells \f$(\gamma_1, \gamma2)\f$ of dimension q / q+1, valid for A (ie.\ such that \f$\langle \partial_{q+1}(\gamma_2), \gamma_1 \rangle\f$ invertible).
+ * It returns a vector of such pairs.
+ *
+ * \param q Lower dimension of the pair.
+ * \param found Reference to a %Boolean variable. The method sets `found` to `true` if a valid pair is found, `false` otherwise.
+ */
+virtual std::vector<Cell_pair> find_pairs_A(int q, bool &found) const;
+
+/**
+ * \brief Finds *all* valid Cell_pair for A containing `gamma` (a cell of dimension `q`)
+ *
+ * The function searches all critical cells \f$\gamma'\f$ such that one of the following conditions holds:
+ * - \f$\gamma'\f$ has dimension q+1 and \f$(\gamma, \gamma')\f$ is valid for A (ie.\ such that \f$\langle \partial_{q+1}(\gamma'), \gamma \rangle\f$ invertible),
+ * - \f$\gamma'\f$ has dimension q-1 and \f$(\gamma', \gamma)\f$ is valid for A (ie.\ such that \f$\langle \partial_{q}(\gamma), \gamma' \rangle\f$ invertible).
+ * It returns a vector of such pairs.
+ *
+ * \param q Dimension of the cell `gamma`.
+ * \param found Reference to a %Boolean variable. The method sets `found` to `true` if a valid pair is found, `false` otherwise.
+ * \param gamma Index of a cell to pair.
+ */
+virtual std::vector<Cell_pair> find_pairs_A(int q, bool &found, size_t gamma) const;
+
+
+/// @}
+
+
+/// \name HDVF functions and operations
+/// @{
+
+/*!
+ * \brief A operation: pairs two critical cells.
+ *
+ * A pair of critical cells \f$(\gamma_1, \gamma_2)\f$ of respective dimension q and q+1 is valid for A if \f$\langle \partial_{q+1}(\gamma_2), \gamma_1 \rangle\f$ is invertible. After the `A()` operation, \f$\gamma_1\f$ becomes PRIMARY, \f$\gamma_2\f$ becomes SECONDARY. The A functions updates the reduction accordingly (in time \f$\mathscr O(n^2)\f$).
+ */
+void A(size_t gamma1, size_t gamma2, int q);
+
+/*!
+ * \brief Computes a perfect HDVF.
+ *
+ * As long as valid pairs for A exist, the function selects the first available pair (returned by `find_pair_A()`) and applies the corresponding `A()` operation.
+ * If the `IntegralDomainWithoutDivision` of coefficients is a field, this operation always produces a perfect HDVF (ie.\ the reduced boundary is null and the reduction provides homology and cohomology information).
+ *
+ * If the HDVF is initially not trivial (some cells have already been paired), the function completes it into a perfect HDVF.
+ *
+ * \param verbose If `true`, all intermediate reductions are printed out.
+ *
+ */
+std::vector<Cell_pair> compute_perfect_hdvf(bool verbose = false);
+
+/*!
+ * \brief Computes a random perfect HDVF.
+ *
+ *As long as valid pairs for A exist, the function selects a random pair (among pairs returned by `find_pairs_A()`) and applies the corresponding `A()` operation.
+ * If the `IntegralDomainWithoutDivision` of coefficients is a field, this operation always produces a perfect HDVF (ie.\ the reduced boundary is null and the reduction provides homology and cohomology information).
+ *
+ * If the HDVF is initially not trivial (some cells have already been paired), the function randomly completes it into a perfect HDVF.
+ *
+ * \param verbose If `true`, all intermediate reductions are printed out.
+ */
+std::vector<Cell_pair> compute_rand_perfect_hdvf(bool verbose = false);
+/// @}
+
+/// \name Getters
+/// @{
+
+/*!
+ * \brief Gets cells with a given `PSC_flag` in any dimension.
+ *
+ * The function returns in each dimension the vector of cells with a given `PSC_flag`.
+ */
+std::vector<std::vector<size_t> > flag (PSC_flag flag) const;
+
+/*!
+ * \brief Gets cells with a given `PSC_flag` in dimension `q`.
+ *
+ * The function returns the vector of cells of dimension `q` with a given `PSC_flag`.
+ */
+std::vector<size_t> psc_flags (PSC_flag flag, int q) const;
+
+/*!
+ * \brief Gets the `PSC_flag` of the cell `tau` in dimension `q`.
+ */
+PSC_flag psc_flag (int q, size_t tau) const;
+
+/*!
+ * \brief Gets HDVF option.
+ */
+int hdvf_opts ();
+
+/*!
+ * \brief Gets the row-major matrix of \f$f\f$ (from the reduction associated to the HDVF).
+ */
+const Row_matrix& matrix_f (int q) const;
+
+/*!
+ * \brief Gets the column-major matrix of \f$g\f$ (from the reduction associated to the HDVF).
+ */
+const Column_matrix& matrix_g (int q) const;
+
+/*!
+ * \brief Gets the column-major matrix of \f$h\f$ (from the reduction associated to the HDVF).
+ */
+const Column_matrix& matrix_h (int q) const;
+
+/*!
+ * \brief Gets the column-major matrix of \f$\partial'\f$, reduced boundary operator (from the reduction associated to the HDVF).
+ */
+const Column_matrix& matrix_dd (int q) const;
+
+/*!
+ * \brief Tests if a HDVF is perfect.
+ *
+ * The function returns `true` if the reduced boundary matrix is null and `false` otherwise.
+ */
+    bool is_perfect_hdvf();
+
+/// @}
+
+/// \name Output functions
+/// @{
+
+/*!
+ * \brief Writes the matrices of the reduction.
+ *
+ * Writes the matrices of the reduction (that is \f$f\f$, \f$g\f$, \f$h\f$, \f$\partial'\f$ the reduced boundary).
+ *
+ * By default, writes the complex to `std::cout`.
+*/
+std::ostream& write_matrices(std::ostream &out = std::cout) const;
+
+/*! \brief Writes the homology and cohomology reduction information.
+ *
+ * Writes the homology and cohomology reduction information (that is \f$f^*\f$, \f$g\f$ \f$\partial'\f$ the reduced boundary over each critical cell).
+ *
+ * By default, writes the complex to `std::cout`.
+*/
+std::ostream& write_reduction(std::ostream &out = std::cout) const;
+
+/*!
+ *\brief Writes a HDVF and its reduction to a stream.
+ *
+ * Write a HDVF together with the associated reduction (f, g, h, d matrices) to the stream `out` in  `.hdvf` file format.
+ */
+
+std::ostream& write_hdvf_reduction(std::ostream& out) ;
+
+/*!
+ * \brief Reads a HDVF and its reduction from a stream.
+ *
+ * Reads a HDVF and its reduction from  the stream `in` in  `.hdvf` file format, a simple text file format (see for a specification).
+ */
+std::istream& read_hdvf_reduction(std::istream& in) ;
+
+/*!
+ * \brief Exports primary/secondary/critical labels (in particular for vtk export).
+ *
+ * The method exports the labels of every cell in each dimension.
+ *
+ * \return A vector containing, for each dimension, the vector of labels by cell index.
+ */
+virtual std::vector<std::vector<int> > psc_labels () const;
+
+/*!
+ * \brief Returns a chain containing the homology generator associated to `cell_index` (critical cell) of dimension  `q` (in particular for vtk export).
+ *
+ * The method exports the chain \f$g(\sigma)\f$ for \f$\sigma\f$ the cell of index `cell_index` and dimension `q`.
+ *
+ * \return A column-major chain.
+ */
+virtual Column_chain homology_chain (size_t cell_index, int q) const;
+
+/*!
+ * \brief Returns a chain containing the cohomology generators associated to `cell_index` (critical cell) of dimension  `q` (in particular for vtk export).
+ *
+ * The method exports:
+ * - the chain \f$f^\star(\sigma)\f$ for \f$\sigma\f$ the cell of index `cell_index` and dimension `q`,
+ * - or the co-faces of this chain if `co_faces` is true (sometimes easier to view cohomology generators)
+ *
+ * Below, a sample mesh with, (left) homology generators, (right) two examples of cohomology generators (corresponding generators/co-generators bear similar colours):
+ *
+ * <img src="HDVF_dtorus_homs.png" align="center" width=25%/>
+ * <img src="HDVF_dtorus_cohom1.png" align="center" width=25%/>
+ * <img src="HDVF_dtorus_cohom2.png" align="center" width=25%/>
+ *
+ * The same generators displayed through their co-faces:
+ *
+ * <img src="HDVF_dtorus_cohom1_co.png" align="center" width=25%/>
+ * <img src="HDVF_dtorus_cohom2_co.png" align="center" width=25%/>
+ *
+ * All homology / cohomology generators:
+ *
+ *<img src="HDVF_dtorus_all.png" align="center" width=30%/>
+ *
+ * \return A column-major chain.
+ */
+virtual Column_chain cohomology_chain (size_t cell_index, int dim, bool co_faces = false) const;
+
+/// @}
+
+};

HDVF/doc/HDVF/Concepts/HDVFTraits.h

@@ -0,0 +1,29 @@
+/*!
+\ingroup PkgHDVFConcepts
+\cgalConcept
+
+The concept `HDVFTraits` describes the requirements for geometric traits classes used in this package.
+
+\cgalHasModelsBegin
+\cgalHasModelsBare{`CGAL::Homological_discrete_vector_field::Hdvf_traits_2<K>`}
+\cgalHasModelsBare{`CGAL::Homological_discrete_vector_field::Hdvf_traits_3<K>`}
+\cgalHasModelsBare{`CGAL::Homological_discrete_vector_field::Hdvf_traits_d<K>`}
+\cgalHasModelsEnd
+
+@todo add requirements on Point, Vector, FT, Bbox
+*/
+class HDVFTraits
+{
+public:
+    /* \brief Type tag encoding the dimension of the data. Must be a `CGAL::Dimension_tag<>` */
+    typedef unspecified_type Dimension ;
+
+    /* \brief Geometric kernel type. Must be a model of the `Kernel` concept. */
+    typedef unspecified_type Kernel ;
+    typedef unspecified_type Point ;
+    typedef unspecified_type Vector ;
+    typedef typename Kernel::FT FT ;
+
+    /* \brief Type of bounding boxes in the data dimension. */
+    typedef unspecified_type Bbox ;
+};

HDVF/doc/HDVF/Concepts/SparseChain.h

@@ -0,0 +1,311 @@
+/*!
+\ingroup PkgHDVFConcepts
+\cgalConcept
+
+The concept `SparseChain` describes the requirements for sparse vectors (called *sparse chains* in homology) optimized for topological computations.
+More precisely, `SparseChain` provides all the operations on chains required by the `SparseMatrix` concept.
+
+ `SparseChain` encodes non zero coefficients of a sparse vector.
+
+ `SparseChain`  can be either row or column vectors. The following constants, called `StorageFormat`, encode the direction of sparse chains (and sparse matrices).
+ - `CGAL::OSM::COLUMN` for column-major chains and matrices (which is the default),
+ - `CGAL::OSM::ROW` for row-major chains and matrices.
+
+
+\cgalHasModelsBegin
+\cgalHasModelsBare{`CGAL::OSM::Sparse_chain<CoefficientRing, StorageFormat>`}
+\cgalHasModelsEnd
+
+ \sa `IntegralDomainWithoutDivision`
+ \sa `SparseMatrix`
+*/
+
+class SparseChain {
+
+public:
+    /// \name Types
+    /// @{
+
+    /*!
+     * \brief Type of coefficients stored in the matrix (a model of `IntegralDomainWithoutDivision`).
+     */
+    typedef unspecified_type Coefficient_ring;
+
+    /*!
+     * \brief Matrix and chain storage format (either CGAL::OSM::ROW or CGAL::OSM::COLUMN).
+     */
+    typedef int Storage_format;
+
+    /*!
+     * \brief Sparse chain iterator type.
+     */
+    typedef unspecified_type iterator;
+
+    /*!
+     * \brief Sparse chain constant iterator type.
+     */
+    typedef unspecified_type const_iterator;
+
+    /// @}
+
+    /// \name Creation, filling
+    /// @{
+
+    /*!
+     * \brief Creates new empty sparse chain.
+     *
+     * Creates a sparse chain encoding an empty linear combination of cells.
+     */
+    SparseChain();
+
+    /*!
+     * \brief Creates new empty sparse chain (ie.\ zero-chain) of given size.
+     *
+     * Constructor with size, initializes an empty sparse chain encoding a linear combination of cells with all coefficients null.
+     */
+    SparseChain(size_t chainSize);
+
+    /*!
+     * \brief Creates new sparse chain by copy.
+     *
+     * Copy constructor, initialize a sparse chain from an existing sparse chain of same storage format.
+     */
+    SparseChain(const SparseChain &other);
+
+    /*!
+     * \brief Assigns to other chain.
+     *
+     * Assign to other chain coefficient-wise, equivalent to copying it.
+     *
+     * SparseChain must have the same `Coefficient_ring`.
+     */
+    SparseChain& operator=(const SparseChain &other);
+
+    /// @}
+
+    /// \name Matrix informations and iterators
+    /// @{
+
+    /*!
+     * \brief Returns the dimension of the basis (that is, size of the chain).
+     */
+    size_t dimension() const;
+
+    /*!
+     * \brief Iterator to the beginning of the chain.
+     *
+     * The function returns an iterator to the first non zero index.
+     */
+    iterator begin() noexcept;
+
+    /*!
+     * \brief Constant iterator to the beginning of the chain.
+     *
+     * The function returns a constant iterator to the first non zero index.
+     */
+    const_iterator begin() const noexcept;
+
+    /*!
+     * \brief Constant iterator to the beginning of the chain.
+     *
+     * The function returns a constant iterator to the first non zero index.
+     */
+    const_iterator cbegin() const noexcept;
+
+    /*!
+     * \brief Past-the-end iterator of the chain.
+     *
+     * The function returns an iterator past the end of the chain.
+     */
+    iterator end() noexcept;
+
+    /*!
+     * \brief Constant past-the-end iterator of the chain.
+     *
+     * The function returns a constant iterator past the end of the chain.
+     */
+    const_iterator end() const noexcept;
+
+    /*!
+     * \brief Constant past-the-end iterator of the chain.
+     *
+     * The function returns a constant iterator past the endi of the chain.
+     */
+    const_iterator cend() const noexcept;
+
+    /// @}
+
+    /// \name Output
+    /// @{
+
+    /*!
+     * \brief Inserts `chain` in the output stream.
+     */
+    friend std::ostream& operator<<(std::ostream &stream, const SparseChain& chain);
+
+    /// @}
+
+    /// \name Linear algebra operators
+    /// @{
+
+    /*!
+     * \brief Adds two chains together.
+     *
+     * Add two chains and return the result in a new matrix.
+     * Chains must have the same `Coefficient_ring` and the same storage format.
+     */
+    friend SparseChain operator+(const SparseChain &first, const SparseChain &second);
+
+    /*!
+     * \brief Subtracts a chain from another chain.
+     *
+     * Subtract two chains and return the result in a new matrix.
+     * Chains must have the same `Coefficient_ring` and the same storage format.
+     */
+    friend SparseChain operator-(const SparseChain &first, const SparseChain &second);
+
+    /*!
+     * \brief Applies multiplication on each coefficient.
+     */
+    friend SparseChain operator*(const Coefficient_ring& lambda, const SparseChain &chain);
+
+    /*!
+     * \brief Applies multiplication on each coefficient.
+     */
+    friend SparseChain operator*(const SparseChain &_chain, const Coefficient_ring& lambda);
+
+    /*!
+     * \brief Performs matrix multiplication between a column chain and a row chain and returns a column-based matrix.
+     *
+     * Generate a column-based matrix from the matrix multiplication and return it.
+     * \pre chains must have the same `Coefficient_ring`.
+     * \pre `column.is_column()`  and `row.is_row()`  must be `true`
+     */
+    friend SparseMatrix operator*(const SparseChain &column, const SparseChain &row);
+
+    /*!
+     * \brief Performs matrix multiplication between a column chain and a row chain and returns a row-based matrix.
+     *
+     * Generate a row-based matrix from the matrix multiplication and return it.
+     * \pre chains must have the same `Coefficient_ring`.
+     * \pre `column.is_column()`  and `row.is_row()`  must be `true`
+     */
+    friend SparseMatrix operator%(const Sparse_chain &column, const Sparse_chain &row);
+
+    /*!
+     * \brief Performs dot product between a row chain and a column chain.
+     *
+     * \pre chains must have the same `Coefficient_ring`.
+     * \pre `row.is_row()` and `column.is_column()` must be `true`
+     */
+    friend Coefficient_ring operator*(const Sparse_chain &row, const Sparse_chain &column);
+
+    /*!
+     * \brief Adds a chain to `this`.
+     *
+     * Add a chain to `this`.
+     * Chains must have the same `Coefficient_ring` and the same storage format.
+     */
+    SparseChain& operator+=(const SparseChain &_other);
+
+    /*!
+     * \brief Subtracts a chain from `this`.
+     *
+     * Subtract a chain from `this`.
+     * Chains must have the same `Coefficient_ring` and the same storage format.
+     */
+    SparseChain& operator-=(const SparseChain &_other);
+
+    /*!
+     * \brief Applies multiplication on each coefficient of `this`.
+     */
+    SparseChain& operator*=(const Coefficient_ring& lambda);
+
+    /*!
+     * \brief Transposes a sparse chain.
+     *
+     * The result is a chain with storage format switched between column and row.
+     */
+    SparseChain transpose();
+
+    /// @}
+
+    /// \name Access and modifications
+    /// @{
+
+    /*!
+     * \brief Compares two chains.
+     */
+    bool operator==(const SparseChain& other_chain);
+
+    /*!
+     * \brief Gets the value of a coefficient of the chain.
+     */
+    Coefficient_ring operator[](size_t index);
+
+    /*!
+     * \brief Gets the value of a coefficient of the chain.
+     */
+    Coefficient_ring get_coefficient(size_t index) const ;
+
+    /**
+     * \brief Sets a given coefficient of the chain.
+     */
+    void set_coefficient(size_t index, Coefficient_ring d);
+
+    /*!
+     * \brief Checks if a coefficient of the chain is null.
+     */
+    bool is_null(size_t index) const;
+
+    /*!
+     * \brief Checks if the chain is null.
+     */
+    bool is_null() const;
+
+    /*!
+     * \brief Gets a sub-chain from the chain.
+     *
+     * Return a new chain where all coefficients of indices provided in the vector are removed.
+     */
+    friend SparseChain operator/(const SparseChain &chain, const std::vector<size_t> &indices);
+
+    /*!
+     * \brief Gets a sub-chain from the chain.
+     *
+     * Return a new chain where the coefficients at a given index is removed.
+     */
+    friend SparseChain operator/(const SparseChain &chain, size_t indices);
+
+    /*!
+     * \brief Restricts the chain to a sub-chain by removing indices.
+     *
+     * Removes all indices provided in the vector from the chain. Return a reference to the modified chain.
+     */
+    SparseChain& operator/=(const std::vector<size_t> &indices);
+
+    /**
+     * \brief Restricts the chain to a sub-chain by removing a given index.
+     *
+     * Removes the index provided from the chain. Return a reference to the modified chain.
+     */
+    SparseChain& operator/=(size_t index);
+
+    /**
+     * \brief Removes all coefficients from the chain.
+     *
+     * The function comes to set all coefficients to zero.
+     */
+    void nullify();
+
+    /*!
+     * \brief Checks if chain is a column.
+     */
+    bool is_column() const;
+
+    /*!
+     * \brief Checks if chain is a row.
+     */
+    bool is_row() const;
+    /// @}
+};

HDVF/doc/HDVF/Concepts/SparseMatrix.h

@@ -0,0 +1,416 @@
+/*!
+\ingroup PkgHDVFConcepts
+\cgalConcept
+
+The concept `SparseMatrix` describes the requirements for sparse matrices optimized for topological computations.
+Traditionally, sparse matrices data structures encode non zero coefficients of (sparse) matrices in order to optimize either matrices memory footprint, or linear algebra operations (which usually comes to optimize iterators over non zero coefficients and access to coefficients). However, topological operations require slightly different features:
+ - fast access to row or columns of matrices (which are actually the images under the application encoded by the matrix)
+ - fast block operations (especially along row or columns)
+
+ The `SparseMatrix` concept describes requirements for such sparse matrix. It relies on the model of `SparseChain` which encodes sparse row or column vectors. Matrices are either column major or row major (hence they either store column sparse chains or row sparse chains).
+
+ The following constants, called `StorageFormat`, encode the major direction of both sparse chains and sparse matrices.
+ - `CGAL::OSM::COLUMN` for column-major chains and matrices (which is the default),
+ - `CGAL::OSM::ROW` for row-major chains and matrices.
+
+ For instance, given the \f$5\times 4\f$ matrix:
+ \f[
+ A = \left(\begin{array}{cccc}
+ 1 & \cdot & \cdot & \cdot \\
+ -1 & \cdot & 2 & \cdot\\
+ \cdot & \cdot & 1 & \cdot \\
+ \cdot & \cdot & \cdot & \cdot \\
+ \cdot & \cdot & \cdot & \cdot \\
+ \end{array}\right)
+ \f]
+ where \f$\cdot\f$ means \f$0\f$.
+ - A column-major representation of \f$A\f$ is: \f$[0\mapsto c_0, 2\mapsto c_2]\f$ with the column-chains: \f$c_0 = [0\mapsto 1, 1\mapsto -1]\f$ and \f$c_2 = [1\mapsto 2, 2\mapsto 1]\f$.
+ - A row-major representation of \f$A\f$ is: \f$[0\mapsto c_0, 1\mapsto c_1, 2\mapsto c_2]\f$ with the row-chains: \f$c_0 = [0\mapsto 1]\f$ and \f$c_1 = [0\mapsto -1, 2\mapsto 2]\f$ and \f$c_2 = [2\mapsto 1]\f$.
+
+\cgalHasModelsBegin
+\cgalHasModelsBare{`CGAL::OSM::Sparse_matrix<IntegralDomainWithoutDivision, StorageFormat>`}
+\cgalHasModelsEnd
+
+ \sa `IntegralDomainWithoutDivision`
+ \sa `SparseChain`
+*/
+
+class SparseMatrix {
+
+public:
+    /// \name Types
+    /// @{
+
+    /*!
+     * \brief Type of coefficients stored in the matrix (a model of `IntegralDomainWithoutDivision`).
+     */
+    typedef unspecified_type Coefficient_ring;
+
+    /*!
+     * \brief Matrix and chain storage format (either `ROW` or `COLUMN`).
+     */
+    typedef int Storage_format;
+
+    /*!
+     * \brief A data structure storing indices of non empty chains.
+     */
+    typedef unspecified_type Non_zero_chain_indices;
+
+    /*!
+     * \brief Type of the chains stored in the matrix.
+     */
+    typedef unspecified_type Matrix_chain;
+
+    /// @}
+
+
+    /// \name Creation, Filling
+    /// @{
+
+    /*!
+     * \brief Creates an empty new sparse matrix object.
+     *
+     * Default constructor, initialize an empty matrix of type `StorageFormat` with coefficients of type `Coefficient_ring`.
+     * The default matrix size is 0x0.
+     */
+    SparseMatrix() ;
+
+    /*!
+     * \brief Creates a new sparse matrix object with given rows/columns sizes.
+     *
+     * Constructor with sizes, initialize an empty matrix of type `StorageFormat` with coefficients of type `Coefficient_ring` and a given size along rows/columns.
+     */
+    SparseMatrix(const size_t rowCount, const size_t columnCount) ;
+
+    /**
+     * \brief Creates a new sparse matrix from another sparse matrix object (with possibly a different `StorageFormat`).
+     *
+     * Copy constructor, initialize a sparse matrix of same sizes, containing the same coefficients (but not necessarly of the same `StorageFormat`).
+     * If types are different, the constructor performs conversion.
+     */
+    SparseMatrix(const SparseMatrix& otherToCopy);
+
+    /*!
+     * \brief Assigns to other matrix.
+     *
+     * Assign to other matrix coefficient-wise, equivalent to copying it.
+     *
+     * Matrices must have the same type.
+     */
+    SparseMatrix& operator=(const SparseMatrix& _otherToCopy);
+
+    /**
+     * \brief Cleans a sparse matrix (set all coefficients to zero).
+     *
+     * Empty all structures of the sparse matrix.
+     *
+     */
+    void nullify();
+
+    /// @}
+
+    /// \name Matrix informations and iterators
+    /// @{
+
+    /**
+     * \brief Tests if a sparse matrix is null.
+     *
+     * The function returns `true` if the sparse matrix is null (that is, empty) and `false` otherwise.
+    */
+    bool is_null();
+
+    /**
+     * \brief Gets the matrix size.
+     *
+     * The matrix size as a row/column pair.
+     */
+    std::pair<size_t, size_t> dimensions() const;
+
+    /**
+     * \brief Iterator to the beginning of the chain indices (visited by increasing indices).
+     *
+     * The function returns an iterator to the beginning of the (non zero) chain indices.
+     */
+    inline Non_zero_chain_indices::iterator begin() const noexcept;
+
+    /**
+     * \brief Past-the-end iterator of the chain indices (visited by increasing indices).
+     *
+     * The function returns an iterator past the end of the chain indices.
+     */
+    inline Non_zero_chain_indices::iterator end() const noexcept;
+
+    /**
+     * \brief Reverse iterator to the beginning of the chain indices (visited by decreasing indices).
+     *
+     * The function returns a reverse iterator to the beginning of the (non zero) chain indices.
+     */
+    inline Non_zero_chain_indices::iterator reverse_begin() const noexcept;
+
+    /**
+     * \brief Past-the-end reverse iterator of the chain indices (visited by decreasing indices).
+     *
+     * The function returns a past-the-end reverse iterator of the chain indices.
+     */
+    inline Non_zero_chain_indices::iterator reverse_end() const noexcept;
+
+    /// @}
+
+    /// \name Output
+    /// @{
+
+    /**
+     * \brief Inserts `matrix` in the output stream.
+     */
+    friend std::ostream& operator<<(std::ostream &_stream, const SparseMatrix &matrix);
+
+    /**
+     * \brief Inserts `matrix` in an output stream.
+     */
+    friend std::ostream& write_matrix (const SparseMatrix& matrix, std::ostream& out);
+
+    /**
+     * \brief Extracts a sparse matrix from an input stream.
+     */
+
+    template <typename _CT>
+    friend std::istream& read_matrix (SparseMatrix& matrix, std::istream& in);
+    /// @}
+
+    /// \name Linear algebra operators
+    /// @{
+
+    /**
+     * \brief Adds a matrix and assign.
+     *
+     * Adds each coefficient of the matrix together and stores the result in `this`.
+     * Matrices must have the same `Coefficient_ring` but can have different `StorageFormat`.
+     */
+    friend SparseMatrix& operator+=(SparseMatrix &matrix, const SparseMatrix &other);
+
+    /**
+     * \brief Adds two matrices together.
+     *
+     * Adds each coefficient of the matrices together and returns a new matrix (of the same type as `first`) representing the result (when possible, prefer `+=` for efficiency).
+     * Matrices must have the same `Coefficient_ring` but can have different `StorageFormat`.
+     */
+    friend SparseMatrix operator+(const SparseMatrix &first, const SparseMatrix &second);
+
+
+    /**
+     * \brief Subtracts a matrix and assign.
+     *
+     * Subtracts each coefficient of the matrix together and stores the result in `matrix`.
+     * Matrices must have the same `Coefficient_ring` but can have different `StorageFormat`.
+     */
+    SparseMatrix& operator-=(SparseMatrix &matrix, const SparseMatrix &other);
+
+    /**
+     * \brief Subtracts two matrices together.
+     *
+     * Subtracts each coefficient of the matrix together and returns a new matrix (of the same type as `first`) representing the result (when possible, prefer `-=` for efficiency).
+     * Matrices must have the same `Coefficient_ring` but can have different `StorageFormat`.
+     */
+    friend SparseMatrix operator-(const SparseMatrix &first, const SparseMatrix &second);
+
+    /**
+     * \brief Computes the negative of a matrix (unary operator).
+     *
+     * \return The resulting matrix.
+     */
+    friend SparseMatrix operator-(const SparseMatrix& matrix);
+
+    /**
+     * \brief Applies factor on each coefficients into a new matrix.
+     *
+     * This method creates a new matrix obtained by multiplying the matrix by a scalar factor `lambda`. If `lambda` is zero, the function comes to nullify the matrix (when possible, prefer `*=` for efficiency).
+     */
+    friend SparseMatrix operator*(const Coefficient_ring& lambda, const SparseMatrix &matrix);
+
+    /**
+     * \brief Applies factor on each coefficients into a new matrix.
+     *
+     * This method creates a new matrix obtained by multiplying the matrix by a scalar factor `lambda`. If `lambda` is zero, the function comes to nullify the matrix (when possible, prefer `*=` for efficiency).
+     */
+    friend SparseMatrix operator*(const SparseMatrix &matrix, const Coefficient_ring& lambda);
+
+    /**
+     * \brief Applies factor on each coefficient and assign.
+     *
+     * This method multiplies the matrix by a scalar factor `lambda`.
+     */
+    SparseMatrix& operator*=(const Coefficient_ring& lambda);
+
+    /**
+     * \brief Multiplies a matrix and assign.
+     *
+     * Multiply each coefficient of the matrix together and stores the result in `matrix`.
+     * Matrices must have the same `Coefficient_ring` but can have different `StorageFormat`.
+     */
+    friend SparseMatrix& operator*=(SparseMatrix &matrix, const SparseMatrix &other);
+
+    /**
+     * \brief Performs multiplication between matrices and returns a new column-major matrix.
+     *
+     * Perform standard linear algebra product between matrices and returns a new column-major matrix (when possible, prefer `*=` for efficiency). Matrices must have the same `Coefficient_ring` but can have different `StorageFormat`. Efficiency of the product depends of `StorageFormat` (when possible, prefer row-major by column-major products).
+     */
+    friend SparseMatrix operator*(const SparseMatrix &first, const SparseMatrix &second);
+
+    /**
+     * \brief Performs multiplication between matrices and returns a new row-major matrix.
+     *
+     * Perform standard linear algebra product between matrices and returns a new row-major matrix (when possible, prefer `*=` for efficiency). Matrices must have the same `Coefficient_ring` but can have different `StorageFormat`. Efficiency of the product depends of `StorageFormat`.
+     */
+    friend SparseMatrix operator%(const SparseMatrix &first, const SparseMatrix &second);
+
+    /**
+     * \brief Performs multiplication between a matrix and a column-based chain.
+     *
+     * Generates a column-based chain from the matrix multiplication and returns it.
+     *
+     * \pre the matrix and the chain must have the same `Coefficient_ring`.
+     * \pre `column.is_column()` must be `true`
+     */
+    friend SparseChain operator*(const SparseMatrix& matrix, const SparseChain& column);
+
+    /**
+     * \brief Performs multiplication between a row-based chain and a matrix.
+     *
+     * Generates a row-based chain from the matrix multiplication and returns it.
+     *
+     * \pre the matrix and the chain must have the same `Coefficient_ring`.
+     * \pre `row.is_row()` must be `true`
+     */
+    friend SparseChain operator*(const SparseChain& row, const SparseMatrix& matrix);
+
+    /**
+     * \brief Transposes a matrix.
+     *
+     * The function returns a new matrix where the `StorageFormat` is changed.
+     */
+    SparseMatrix transpose();
+
+
+    /// @}
+
+    /// \name Access and blocks operations
+    /// @{
+
+    /**
+     * \brief Gets a chain from a const matrix.
+     *
+     * If the matrix is column-major, returns the `i`th column, and if the matrix is row-major, returns the `i`th row.
+     */
+    Matrix_chain operator[](size_t index) const;
+
+    /**
+     * \brief Sets a given coefficient.
+     *
+     * Assign the scalar `d` to the coefficient on row `i` and column `j`.
+     */
+    friend void set_coefficient(SparseMatrix& matrix, size_t i, size_t j, const Coefficient_ring d);
+
+    /**
+     * \brief Gets a given coefficient.
+     *
+     * Returns the coefficient on row `i` and column `j` of the matrix.
+     */
+    friend Coefficient_ring get_coefficient(const SparseMatrix& matrix, size_t i, size_t j);
+
+    /**
+     * \brief Gets the value of the column at a given `index` from the matrix (whatever the `StorageFormat` of the matrix).
+     *
+     * For column-matrices, it is equivalent to `operator[]`, for row-matrices a traversal of the matrix is required (in \f$\mathcal O(n)\f$).
+     *
+     *\returns a column chain.
+     */
+    friend SparseChain get_column(const SparseMatrix &matrix, size_t index);
+
+    /**
+     * \brief Gets the value of the row at a given `index` from the matrix (whatever the `StorageFormat` of the matrix).
+     *
+     * For row-matrices, it is equivalent to `operator[]`, for column-matrices a traversal of the matrix is required (in \f$\mathcal O(n)\f$).
+     *
+     *\returns a row chain.
+     */
+    friend SparseChain get_row(const SparseMatrix &matrix, size_t index);
+
+    /**
+     * \brief Gets a constant reference over the column of  index`i` from a column matrix.
+     * \pre `matrix.is_column()` must be `true`.
+     *
+     * \returns a constant reference to a column chain.
+     */
+    const SparseChain& cget_column(const SparseMatrix& matrix, size_t i);
+
+    /**
+     * \brief Gets a constant reference over the row of  index`i` from a row matrix.
+     * \pre `matrix.is_row()` must be `true`.
+     *
+     * \returns a constant reference to a row chain.
+     */
+    const SparseChain& cget_row(const SparseMatrix& matrix, size_t i);
+
+    /**
+     * \brief Sets a column in the matrix (whatever the `StorageFormat` of the matrix).
+     *
+     * Set the `i`th column of `matrix` to `column`.
+     * For column-matrices, it should be equivalent to an assignment, however, for row-matrices, a traversal of the matrix is required (in \f$\mathcal O(n)\f$).
+     * \pre `column.is_column()` must be `true`
+     */
+    void set_column(SparseMatrix &matrix, size_t i, const SparseChain &column);
+
+    /**
+     * \brief Sets a row in `matrix` (whatever the `StorageFormat` of the matrix).
+     *
+     * Set the `i`th row of `matrix` to `chain`.
+     * For row-matrices, it should be equivalent to an assignment, however, for column-matrices, a traversal of the matrix is required (in \f$\mathcal O(n)\f$).
+     * \pre `row.is_row()` must be `true`
+     */
+    void set_row(SparseMatrix &matrix, size_t i, const Sparse_chain &row);
+
+    /**
+     * \brief Gets a submatrix from the matrix.
+     *
+     * Nullifies the chain of index `i` along the major direction of a copy of the matrix amd returns it.
+     */
+    friend SparseMatrix operator/(const SparseMatrix &matrix, size_t i);
+
+    /**
+     * \brief Gets a submatrix from the matrix and assign.
+     *
+     * Removes (along the major dimension) all indices provided in the vector `indices` from the matrix and returns it.
+     */
+    SparseMatrix& operator/=(const std::vector<size_t> &indices);
+
+    /**
+     * \brief Gets a submatrix from the matrix and assign.
+     *
+     * Removes (along the major dimension) the chain of index `i` from the matrix and returns it.
+     */
+    SparseMatrix& operator/=(size_t i);
+
+    /**
+     * \brief Nullifies a column from the matrix.
+     *
+     * Removes column of index `i` whatever the `StorageFormat` of the matrix. For column matrices, it just comes to the `\=` operator and for row matrices, it entails a traversal of the matrix.
+     */
+    friend SparseMatrix& remove_column(SparseMatrix& matrix, size_t index);
+
+    /**
+     * \brief Nullifies a row from the matrix.
+     *
+     * Removes row of index `i` whatever the `StorageFormat` of the matrix. For row matrices, it just comes to the `\=` operator and for column matrices, it entails a traversal of the matrix.
+     */
+    friend SparseMatrix& remove_row(SparseMatrix& matrix, size_t index);
+
+    /**
+     * \brief Nullifies a coefficient of the matrix.
+     *
+     * Removes coefficient on row `i` / column `j` of the matrix.
+     */
+    friend SparseMatrix& remove_coefficient(SparseMatrix& matrix, size_t i, size_t j);
+
+    /// @}
+
+};

HDVF/doc/HDVF/Doxyfile.in

@@ -0,0 +1,3 @@
+@INCLUDE = ${CGAL_DOC_PACKAGE_DEFAULTS}
+
+PROJECT_NAME = "CGAL ${CGAL_DOC_VERSION} - Homological Discrete Vector Fields"

HDVF/doc/HDVF/HDVF.txt

@@ -0,0 +1,790 @@
+namespace CGAL {
+/*!
+
+\mainpage User Manual
+\anchor Chapter_HDVF
+
+\cgalAutoToc
+\author Alexandra Bac
+
+Algebraic topology, introduced more than a century ago, studies topological spaces using algebraic tools derived from them (e.g. groups). Since the 70's and 80's, pioneering works (MacLane, Hilton, Munkres...) paved the way to computational homology, which studies the ``holes'' of an object.
+In particular, Munkres introduced an algorithm computing Betti numbers based on the pseudo-diagonalization of boundary matrices (Smith Normal Form).
+
+In the 90's, persistent homology emerged independently through the work of H. Edelsbrunner et al. (North Carolina), P. Frosini et al. (Bologna), and V. Robins (Colorado), extending the study of homology to sequences of complexes (filtrations).
+This later theory proved to have a wide range of applications in topological data analysis and hence quickly developed and spread in the scientific community.
+In parallel, in the late 80's, R. Forman defined a discrete version of Morse theory for CW-complexes, while F. Sergeraert introduced effective homology theory.
+
+In [AGL, 2017], A. Gonzalez-Lorenzo et al. defined the **homological discrete vector field** (which we refer to as HDVF in the present package), a combinatorial object that generalizes previous approaches.
+Precisely, HDVFs combine discrete Morse and effective homology theory and fully encode homology information (homology groups decomposition, Betti numbers, homology generators) for cellular complexes.
+HDVFs encompass discrete Morse theory, effective homology and Smith Normal form, and: persistent homology can be computed with HDVFs (see [AGL, 2016]), relative homology and Alexander duality (see [AGL, 2025]).
+
+ \cgalFigureBegin{HDVFtwirl,HDVF_twirl.png}
+ HDVF on a twirl model. **Left:** a perfect HDVF over the model; **Middle:** Homology generators associated to this HDVF; **Right:** Corresponding cohomology generators.
+ \cgalFigureEnd
+
+Unlike other main stream state of art approaches and librairies (largely based on persistent homology), HDVFs provide a framework for "geometric" homology and cohomology computation, where holes are not only enumerated, but also "localized" (through homology, cohomology generators, annotations...).
+
+
+
+\section secHDVFdefinitions Overview
+
+The HDVF package is dedicated to homology and cohomology computation for discrete objects and provides tools to compute:
+- basic homological information (such as Betti numbers, ie.\ number of holes per dimension),
+- but also geometric, higher level information (homology and cohomology generators, homology/cohomology annotation, cycles and co-cycles comparison...),
+- moreover HDVFs compute persistent homology (and provide "persistent generators" information) as well as Alexander duality and relative homology.
+
+\section secHDVFHomology Homology/cohomology and related reminders
+
+Homology and cohomology are well know topological invariants which accounts for the holes of an object using commutative groups built from the object. In Mathematics, homology and cohomology can be defined for general Hausdorff topological spaces; we speak then about "singular homology and cohomology". However, such a general definition is not computable and thus, homology/cohomology computation requires a discretization of spaces into cells of various dimensions. Such a discretization is called a **complex**:
+- **simplicial complexes** are intuitively meshes built from vertices, edges, triangles, tetrahedra...
+- **cubical complexes** are built from vertices, edges, squares, cubes...
+- **cellular complexes** are built from vertices, edges, polygons (without holes), polyhedra (without holes)...
+
+Actually, a set of cells is called a complex \f$K\f$ when:
+- it is closed under the "face" operation (ie.\ all the faces of a cell of the complex also belong to the complex)
+- intersections of cells are either empty of they are a proper cell of the complex
+The dimension of \f$K\f$ is the largest dimension of its cells.
+
+ \cgalFigureBegin{HDVFcomplexes,HDVF_cub_simpl_complex.png}
+ Illustration of complexes. **Left:** a simplicial complex of dimension 3; **Right:** A cubical complex of dimension 3.
+ \cgalFigureEnd
+
+ In order to compute homology and cohomology, algebraic objects are built from such "geometric" complexes, namely chain complexes.
+
+\subsection secHDVFChainComplexes Chain complexes
+
+In order to define this algebraic structure, one first need to chose a **ring of coefficients** \f$A\f$  (\f$A\f$ is often taken to be \f$Z\f$ or a \f$\mathbb Z/p\mathbb Z\f$ field).
+
+Given a complex \f$K\f$ of dimension \f$d\f$, a **chain complex** \f$C\f$ is:
+
+\f\[C = \cdots C_{q+1} \xrightarrow{\partial_{q+1}} C_q \xrightarrow{\partial_{q}} C_{q-1} \xrightarrow{\partial_{q-1}} \cdots \to C_1 \xrightarrow{\partial_1} C_0 \xrightarrow{\partial_0 = 0} 0\f\]
+
+
+
+- a sequence of commutative groups, called **chain groups** \f$C_q\f$ for any \f$0\leqslant q \leqslant d\f$ defined as the formal linear combinations of cells with coefficients in \f$A\f$ (or \f$A\f$-module).<br>
+Intuitively, elements of \f$C_q\f$, called **chains** are sums \f$\gamma = \sum_{i=1}^k\lambda_i\cdot \sigma_i\f$ with \f$\lambda_i\in A\f$ and \f$\sigma_i\f$ cells of dimension \f$q\f$, which can be imagined as a set of paths of oriented cells with multiplicity (see \cgalFigureRef{HDVFchain_complexes} where two 1-chains are represented).
+- a sequence of morphisms between chain groups, called **boundary morphisms** \f$\partial_q\,:\, C_q\to C_{q-1}\f$ for any \f$0\leqslant q \leqslant d\f$, "mimicking" algebraically the geometric boundary and satisfying
+
+\f\[\partial_q\partial_{q+1} = 0\f\]
+
+Coming back to the example of \cgalFigureRef{HDVFchain_complexes}, as we will see in section \ref secHDVFsimplicial_chain_complex, the boundary of the edge \f$e_{45}\f$ is the 0-chain \f$v_5-v_4\f$ (ordering is given by vertices indices). Hence, on easily gets that \f$\partial_1(\gamma_1)=0\f$ and \f$\partial_1(\gamma_2)=v_8-v_4\f$, which is a very intuitive result.
+
+ \cgalFigureBegin{HDVFchain_complexes,simplicial_complex_chains.png}
+ Examples of chains representation in a simplicial complex (with \f$A = \mathbb Z\f$). **Pink:** representation the 1-chain (of edges) \f$\gamma_1 = e_{24}+e_{46}-e_{26}\f$; **Blue:** representation of the 1-chain \f$\gamma_2 = e_{45}+e_{56}+e_{67}+e_{78}\f$.
+ \cgalFigureEnd
+
+The definition of the boundary operator depends on the underlying complex, it is detailed in next sections for simplicial and cubical complexes respectively.
+
+ <B>Notation</B>
+
+Given a chain \f$\gamma\in C_q\f$ and \f$\sigma\f$ a \f$q\f$-cell, we denote by \f$\langle\gamma,\sigma\rangle\f$ the coefficient of \f$\sigma\f$ in \f$\gamma\in C_q\f$.
+
+\subsubsection secHDVFsimplicial_chain_complex Simplicial chain complexes
+
+An **abstract simplicial complex** of dimension \f$d\f$ is built from:
+- a set of vertices identified by an index of type `size_t` (cells of dimension 0 or **0-simplices**)
+- for any \f$q\f$ with \f$0< q\leqslant d\f$, cells of dimension \f$q\f$ (or **q-simplices**) are defined by \f$q+1\f$ vertices (denoted \f$\langle v_0,\ldots,v_{q}\rangle\f$ with vertices sorted by *increasing indices*).
+
+A **geometric simplicial complex** equips vertices with coordinates in \f$\mathbb R^q\f$ so that \f$q\f$-simplices are convex hulls of \f$q+1\f$ points.
+
+Given a simplicial complex (abstract or geometric) \f$K\f$, the associated chain complex \f$(C,\partial)\f$, called a **simplicial chain complex** is defined as follows:
+
+- The \f$q\f$-th chain group \f$C_q\f$ is the free \f$A\f$-module generated by \f$q\f$-cells described in the section \ref secHDVFChainComplexes.
+
+- The \f$q\f$-th boundary morphism is defined over the \f$q\f$-simplices (ie.\ the basis) by:
+
+\f\[\partial_q : \begin{array}{rcl}
+C_q & \to & C_{q-1} \\
+\langle v_0,\ldots,v_{q}\rangle & \mapsto & \sum_{i=0}^q (-1)^i\cdot \langle v_0,\ldots,\hat{v_i}, \ldots, v_{q}\rangle \\
+\end{array}
+\f\]
+where \f$\hat{v_i}\f$ means that the \f$i\f$-th vertex has been discarded.
+
+For instance, in our previous example:
+- For the edge \f$\langle 4,5\rangle\f$, we get \f$\partial_1(\langle 4,5\rangle) = \langle 5\rangle - \langle 4\rangle\f$
+- For the triangle \f$\langle 2,4,6\rangle\f$, we get \f$\partial_2(\langle 2,4,6\rangle) = \langle 4,6\rangle - \langle 2,6\rangle + \langle 2,4\rangle\f$
+
+\subsubsection secHDVFcubical_chain_complex Cubical chain complexes
+
+A **cubical complex** of dimension \f$d\f$ is a set of cells defined as cartesian products:
+\f\[\sigma = [n_1,n_1+\delta_1]\times\cdots \times [n_d,n_d+\delta_d]\f\]
+with \f$n_i\in \mathbb N\f$ and \f$\delta_i\in\{0,1\}\f$ for \f$i=1\ldots d\f$.
+The dimension of such a cell is \f$\sum_{i=1}^d \delta_i\f$.
+As stated in previous section, a complex is closed under the "face" operation (if a cell belongs to the complex, so do its faces).
+
+**Khalimsky coordinates** associate to each cell, whatever its dimension, a unique coordinate in \f$\mathbb N^d\f$ defined as:
+\f\[\mathrm{khal}(\sigma) = (2 n_1+\delta_1,\ldots, 2 n_d+\delta_d)\f\]
+The dimension of a cell in this coordinate system is thus the number of odd coordinates.
+
+ \cgalFigureBegin{HDVFcubical_complex_fig,cubical_complex.png}
+ A cubical complex of dimension 3.
+ \cgalFigureEnd
+
+ For instance, in the example above:
+ - \f$v = [2,2] \times [1,1] \times [0,0]\f$ and \f$\mathrm{khal}(\sigma) = (4,2,0)\f$ (cell of dimension 0)
+ - \f$e = [0,1] \times [1,1] \times [0,0]\f$ and \f$\mathrm{khal}(\sigma) = (1,2,0)\f$ (cell of dimension 1)
+ - \f$f = [2,3] \times [0,1] \times [1,1]\f$ and \f$\mathrm{khal}(\sigma) = (5,1,2)\f$ (cell of dimension 2)
+
+Given a cubical complex \f$K\f$, the associated chain complex \f$(C,\partial)\f$, called a **cubical chain complex** is defined as follows:
+- The \f$q\f$-th chain group \f$C_q\f$ is the free \f$A\f$-module generated by \f$q\f$-cells described in the section \ref secHDVFChainComplexes.
+- The \f$q\f$-th boundary morphism is defined over Khalimsky coordinates of cells (ie.\ the basis) as follows:
+\f\[\partial_q(x_1,\ldots,x_d) = \sum_{\substack{i=1\\i\text{ odd}}}^d (-1)^{i+1}\Big((x_1,\ldots, x_i+1,\ldots, x_d) - (x_1,\ldots, x_i-1,\ldots, x_d)\Big)\f]
+
+For instance, in the example above:
+- \f$\partial(e) = (2,2,0)-(0,2,0)\f$
+- \f$\partial(f) = \big((6,1,2)-(4,1,2)\big) - \big((5,2,2)-(5,0,2)\big)\f$
+
+\subsubsection secHDVFcochain Co-chain complex (for cohomology)
+
+In the "mathematical world", singular cohomology is defined on dual groups \f$C^q = C_q^* = \mathrm{Hom}(C_q,A)\f$ with coboundary morphisms: \f$\partial^q\,:\, C^q \to C^{q+1} = \partial_q^*\f$.
+
+However in our discrete setting of finite complexes, a space and its dual can be identified. Hence cohomology can be defined in a more intuitive way:
+
+\f\[C^* = \cdots C^{q+1} \xleftarrow{\partial^{q}} C^q \xleftarrow{\partial^{q-1}} C^{q-1}  \leftarrow \cdots \leftarrow C^2 \xleftarrow{\partial^1} C^1 \xleftarrow{\partial^0} C^0\f\]
+
+- **cochain groups** \f$C^q\f$ for any \f$0\leqslant q \leqslant d\f$ are identified with \f$C_q\f$. Hence chains and cochains are similar.
+- a sequence of morphisms between chain groups, called **coboundary morphisms** \f$\partial^q\,:\,C^q\to C^{q+1}\f$ for any \f$0\leqslant q \leqslant d\f$ with \f$\partial^q = \partial_q^t\f$ (thus, the matrix of the coboundary operator is the transpose of the boundary operator matrix). The coboundary operator "mimicks" algebraically the geometric coboundary (ie.\ cofaces) and satisfies:
+\f\[\partial^q\partial^{q-1} = 0\f\]
+
+The following figure illustrates the coboundary morphism
+
+ \cgalFigureBegin{HDVFcoboundary,simplicial_complex_coboundary.png}
+ Illustration of the coboundary morphism (with \f$A = \mathbb Z\f$) **Pink:** coboundary of vertex \f$v_6\f$ (\f$\partial^0(v_6) = e_{26} + e_{46}+e_{56} - e_{67} - e_{68}\f$); **Blue:** coboundary of edge \f$e_{24}\f$ (\f$\partial^1(e_{24} ) = f_{246} - f_{234}\f$).
+ \cgalFigureEnd
+
+\subsection secHDVFhomology Homology of discrete chain complexes
+
+Homology is a topological invariant characterizing the holes of an object of various dimension. What is a hole?
+- holes of dimension 0 are connected components
+- holes of dimension 1 are tunnels or handles of the object
+- holes of dimension 2 are cavities
+
+ \cgalFigureBegin{HDVFholes_example,holes_example.png}
+ Illustration of the holes of a double torus. **Black:** 1 connected component (dimension 0); **Red:** 4 holes of dimension 1; **Blue** 1 cavity of dimension 2.
+ \cgalFigureEnd
+
+Let us now come back to our simplicial complex example:
+ \cgalFigureBegin{HDVFsimplicial_cycles_example,simplicial_complex_homology.png}
+ Cycles on a simplicial complex. **Blue:** this cycle is the boundary of \f$f_{234}+f_{246}\f$; **Pink:** this cycle is not a boundary.
+ \cgalFigureEnd
+
+This example suggests that holes are delineated by *cycles which are not a boundaries*:
+- **Cycles** of dimension \f$q\f$ are chains in the kernel of \f$\partial_q\f$:
+\f\[Z_q(C) = \ker(\partial_q) \f\]
+- **Boundaries** of dimension \f$q\f$ are chains in the image of \f$\partial_{q+1}\f$
+\f\[B_q(C) = \mathrm{Im}(\partial_{q+1}) \f\]
+- As \f$\partial_{q}\partial_{q+1} = 0\f$, \f$B_q(C)\subseteq Z_q(C)\f$.
+- Hence we can define the quotient \f$Z_q(C)/B_q(C)\f$ which contains exactly "cycles which are not boundaries", ie.\ the \f$q\f$th homology group:
+\f\[H_q(C) = Z_q(C)/B_q(C) \f\]
+
+For finite complexes, these groups are finitely generated and we have the following decomposition (theorem of decomposition of finitely generated groups):
+\f\[H_q(C) = A^{\beta_q}\oplus A/\lambda_1A\oplus \cdots \oplus A/\lambda_k A \f\]
+where:
+- \f$\beta_q\in\mathbb N\f$ is called the **\f$q\f$th Betti number** (giving the number of holes of dimension \f$q\f$);
+- coefficients \f$\lambda_i\f$ are called the **torsion coefficients in dimension \f$q\f$** and satisfy \f$\lambda_i | \lambda_{i+1}\f$.
+
+Let us assume homology groups have no torsion; homology can be computed "at different levels":
+- **level 1:** Computing Betti numbers.
+- **level 2:** Computing \f$\beta_q\f$ cycles generating \f$H_q(C)\f$ called **generators** (hence cycles "delineating" the holes).
+
+Most libraries focus on persistent homology computation for topological data analysis (see section \ref secHDVFpersistent_homology), which is a "level 1" homology computation.
+The **HDVF package is dedicated to a "level 2" homology computation**, which is computationally more complex but provides "geometric" informations about holes (generators, annotations...).
+
+\subsection secHDVFcohomology Cohomology of discrete chain complexes
+
+Starting from a cochain complex:
+
+\f\[C^* = \cdots C^{q+1} \xleftarrow{\partial^{q}} C^q \xleftarrow{\partial^{q-1}} C^{q-1}  \leftarrow \cdots \leftarrow C^2 \xleftarrow{\partial^1} C^1 \xleftarrow{\partial^0} C^0\f\]
+
+one can define:
+- **Cocycles** of dimension \f$q\f$ are (co)chains in the kernel of \f$\partial^q\f$:
+\f\[Z^q(C^*) = \ker(\partial^q) \f\]
+- **Coboundaries** of dimension \f$q\f$ are chains in the image of \f$\partial^{q-1}\f$
+\f\[B^q(C^*) = \mathrm{Im}(\partial^{q-1}) \f\]
+- As \f$\partial^{q}\partial^{q-1} = 0\f$, \f$B^q(C^*)\subseteq Z^q(C^*)\f$.
+- Hence we can define the quotient \f$Z^q(C^*)/B^q(C^*)\f$, ie.\ the \f$q\f$th cohomology group:
+\f\[H^q(C^*) = Z^q(C^*)/B^q(C^*) \f\]
+
+ \cgalFigureBegin{HDVFex_hom_cohom,HDVF_ex_hom_cohom.png}
+ Illustration of homology and cohomology generators for a simple 2D cubical complex \f$C\f$ (computed with HDVFs). **Left:** both generators of \f$H_1(C)\f$; **Middle, right:** both cohomology generators \f$H^1(C)\f$ (for convenience, the dual of cohomology generators has been displayed with bulleted pink lines).
+ \cgalFigureEnd
+
+ \cgalFigureBegin{HDVFex_hom_cohom_cofaces,HDVF_ex_hom_cohom_cofaces.png}
+ Cofaces representation of cohomology generators.
+ \cgalFigureEnd
+
+Figure \cgalFigureRef{HDVFex_hom_cohom} and \cgalFigureRef{HDVFex_hom_cohom_cofaces} illustrate homology and cohomology generators in a simple cubical complex. Actually:
+ - homology generators, intuitively delineate holes;
+ - cohomology generators provide cocycles to open or break the holes (therefore, it is sometimes more convenient to display the cofaces of cohomology generators as in figure \cgalFigureRef{HDVFex_hom_cohom_cofaces}).
+
+Therefore homology and cohomology provide isomorphic (and dual) results. Actually, the following proposition holds in finite dimension:
+\f\[H^q(C^*) \sim H_q(C)\f\]
+
+
+\subsection secHDVFcompHomology Computing homology and cohomology
+
+First results on the computation of discrete homology and cohomology date back to the 80's (see for instance [Munkres]). They are based on a diagonalisation of boundary matrices (called Smith Normal form).
+Nowadays, computational (co)homology approaches can be classified into three branches illustrated in figure \cgalFigureRef{HDVFoverview_comput_hom}.
+
+ \cgalFigureBegin{HDVFoverview_comput_hom,HDVFoverview_comput_hom.png}
+ Overview of approaches in computational homology. HDVFs are in-between discrete Morse theory and effective homology (and they encompass persistent homology, relative homology and tri-partitions).
+ \cgalFigureEnd
+
+\subsubsection secHDVFreduction Effective homology and reductions
+
+Effective homology considers homology computations from a categorical perspective ; its main results are theorems building functors providing a computation of homology for infinite spaces (locally effective spaces) starting from finite ones (effective spaces).
+
+Among the tools of effective homology, **reductions** characterise how to map a (large) chain complex onto a (small) one with isomorphic homology groups. Formally, a reduction from a chain complex \f$C\f$ to a chain complex \f$C'\f$ is a triple of graded maps \f$(f,g,h)\f$:
+\f\[\begin{array}{cccccc}
+\cdots & C_{q+1} & \substack{\xrightarrow{\partial_{q+1}}\\ \xleftarrow[h_{q}]{}} & C_{q} & \substack{\xrightarrow{\partial_{q}}\\ \xleftarrow[h_{q-1}]{}} & \cdots \\
+& f_{q+1}\Big\downarrow \Big\uparrow g_{q+1} & & f_{q}\Big\downarrow \Big\uparrow g_{q} & & \\
+\cdots & C'_{q+1} & \xrightarrow{\partial'_{q+1}} & C'_{q} & \xrightarrow{\partial'_{q}} & \cdots \\
+\end{array}
+\f\]
+such that:
+
+- \f$f\f$ and \f$g\f$ are chain maps (ie.\ commute with \f$\partial\f$ and \f$\partial'\f$)
+- \f$f_q g_q = \mathrm{Id}_{C'_q}\f$ for any \f$q\f$
+- \f$g_q f_q = \mathrm{Id}_{C_q} - \partial_{q+1} h_q - h_q \partial_{q-1}\f$  for any \f$q\f$
+- \f$h_{q+1} h_q = 0\f$, \f$f_{q+1} h_{q} = 0\f$ and \f$h_{q} g_q = 0\f$ for any \f$q\f$
+
+
+\f$C'\f$ is called the **reduced complex** and \f$\partial'\f$ the **reduced boundary**.
+
+Under these assumptions, the homology and cohomology groups of \f$C\f$ and \f$C'\f$ are isomorphic.
+
+A reduction is **perfect** when the reduced boundary is null; then for any \f$q\f$:
+\f\[H_q(C) \sim g(C'_q)\f\]
+in other words, homology groups are isomorphic to the image of the reduced complex by \f$g\f$,
+\f\[H^q(C) \sim f^*(C'_q)\f\]
+and cohomology groups are isomorphic to the image of the reduced complex by \f$f^*\f$.
+
+\subsubsection secHDVFdiscrete_Morse A few words about discrete Morse theory
+
+Discrete Morse theory was introduced by Forman in the 90's as a combinatorial approach for homology computation directly inspired by Morse theory (on compact manifolds).
+He showed that in the discrete context, Morse functions are better defined through their gradient: a vector field defined over the complex connecting cells and one of their co-faces. In this context, gradients of Morse functions can be easily characterized by acyclicity conditions.
+
+More precisely, given a complex \f$K\f$, a **discrete gradient vector field** (called DGVF here after) can be built as follows:
+- consider the Hasse diagram of the complex (connecting cells to their boundary - figure (left))
+
+\subsection secHDVFpersist_hom Persistent homology
+
+\subsection secHDVFrelative_hom Relative homology
+
+\subsection secHDVFAlexander Alexander duality
+
+\section secHDVF The HDVF framework
+
+This section presents a brief introduction to HDVFs, see [AGL, 2016] for more details.
+
+\subsection secHDVFHDVF_definitions Definition of HDVFs
+
+HDVFs are an "in-between" theory:
+- in-between Discrete Morse theory (combinatorial approach) and effective homology (categorical approach)
+- in-between combinatorics and algebra
+but they cannot be reduced to any of these frameworks. As we will see, they provide a bridge in-between them.
+
+Let \f$K\f$ be a complex (simplicial, cubical...) and let \f$(\mathscr C, \partial)\f$ be its associated chain complex.
+
+**Notation** Starting from now on, we will use \f$\mathscr C\f$ for chain complexes and \f$C\f$ for critical cells.
+
+Formally a **HDVF** \f$X=(P,S,C)\f$ is a partition of the cells of \f$K\f$:
+ \f\[K = P\cup S\cup C\f\]
+ with \f$P\f$ primary cells, \f$S\f$ secondary cells and \f$C\f$ critical cells, such that the sub-matrix:
+ \f\[\partial (S)|_P\text{ is invertible.}
+ \f\]
+where \f$\partial (S)|_P\f$ means that the matrix of \f$\partial\f$ is restricted to columns of secondary cells and rows of primary cells.
+
+As stated earlier, this definition is both combinatorial (partition of cells) and algebraic (invertibility of a sub-matrix) but it cannot be simplified or reduced to either of these approaches.
+
+Given such an HDVF \f$X=(P,S,C)\f$, we can define the following associated reduction \f$(f,g,h)\f$:
+
+\f\[\begin{array}{cccccc}
+\cdots & \mathscr C_{q+1} & \substack{\xrightarrow{\partial_{q+1}}\\ \xleftarrow[h_{q}]{}} & \mathscr C_{q} & \substack{\xrightarrow{\partial_{q}}\\ \xleftarrow[h_{q-1}]{}} & \cdots \\
+& f_{q+1}\Big\downarrow \Big\uparrow g_{q+1} & & f_{q}\Big\downarrow \Big\uparrow g_{q} & & \\
+\cdots & C_{q+1} & \xrightarrow{\mathrm d_{q+1}} & C_{q} & \xrightarrow{\mathrm d_{q}} & \cdots \\
+\end{array}
+\f\]
+\f$\mathrm d\f$ is the **reduced boundary** (over critical cells) and the reduction is defined as follows:
+
+- the matrix of \f$h_q\f$ is:
+\f\[\begin{array}{c}
+\ \ \ \ \ \ \ \ \ \ P_q \ \ \ \ S_q \ \ \ \  C_q \\
+\begin{array}{c}
+P_{q+1} \\
+S_{q+1} \\
+C_{q+1} \\
+\end{array}
+\begin{array}{|c|c|c|}
+\hline
+0 & 0 & 0 \\
+\hline
+H_q & 0 & 0 \\
+\hline
+0 & 0 & 0 \\
+\hline
+\end{array}
+\end{array}
+\f\]
+with \f\[H_q = (\partial (S_{q+1})|_{P_q})^{-1}\f\].
+
+- the matrix of \f$f_q\f$ is:
+\f\[\begin{array}{c}
+\ \ \ \ \ \ P_q \ \ \ \ S_q \ \ \ \  C_q \\
+\begin{array}{c}
+C_q \\
+\end{array}
+\begin{array}{|c|c|c|}
+\hline
+F_q & 0 & \mathrm{Id} \\
+\hline
+\end{array}
+\end{array}
+\f\]
+with
+\f\[ F_q = - \partial(S_{q+1})|_{C_q}\cdot H_q\f\]
+
+- the matrix of \f$g_q\f$ is:
+\f\[\begin{array}{c}
+\ \ \ \ \ \ C_q \\
+\begin{array}{c}
+P_q \\
+S_q \\
+C_q \\
+\end{array}
+\begin{array}{|c|}
+\hline
+0 \\
+\hline
+G_q \\
+\hline
+\mathrm{Id} \\
+\hline
+\end{array}
+\end{array}
+\f\]
+with
+\f\[ G_q = - H_{q-1} \cdot \partial(C_{q})|_{P_{q-1}}\f\]
+
+- last, the matrix of the reduced boundary \f$\mathrm d_q\f$ is:
+\f\[\begin{array}{c}
+\ \ \ \ \ \ C_q \\
+\begin{array}{c}
+C_{q-1} \\
+\end{array}
+\begin{array}{|c|}
+\hline
+DD_q \\
+\hline
+\end{array}
+\end{array}
+\f\]
+with
+\f\[ D_q = \partial(C_{q})|_{C_{q-1}} + F_{q} \cdot \partial(C_{q})|_{P_{q-1}} = \partial(C_{q})|_{C_{q-1}} + \partial(S_{q})|_{C_{q-1}} \cdot G_q \f\]
+
+Therefore the HDVF "reduces" the complex to a smaller complex over critical cells.
+
+A HDVF is called **perfect** when the reduced boundary is null, then for any \f$q \in \{0\ldots d\}\f$:
+- \f$\beta_q = | C_q |\f$
+- Homology generators of dimension \f$q\f$ are \f$\{g(\sigma) \,;\, \sigma\in C_q \}\f$
+- Cohomology generators of dimension \f$q\f$ are \f$\{f^*(\sigma) \,;\, \sigma\in C_q \}\f$
+
+Hence a perfect HDVF computes homology and provides a rich homological information.
+- when \f$A\f$ is a field, perfect HDVFs always exist and any HDVF can be completed into a perfect HDVF (see \ref secHDVFcomput)
+- when \f$A\f$ is a ring:
+    - if the torsion part is null, we hypothesise that a perfect HDVF always exist
+    - with torsion, we hypothesise that a torsion-perfect HDVF always exist (where the reduced boundary is limited to torsion)
+Both questions are still an open and we have no counter-examples.
+
+Moreover, the reduction provides various informations on cycles and co-cycles:
+- An **annotation** associates to a cycle its decomposition in the basis of homology generators, up to boundaries: <br>
+Given \f$\gamma\in Z_q(\mathscr C)\f$, \f$f(\gamma)\f$ is an annotation of \f$\gamma\f$:
+\f\[\text{if }f(\gamma) = \sum_i \alpha_i\cdot \sigma_i,\ \ \ \sigma_i\in C_q\f\]
+\f\[\text{then }\ \ \ \gamma = \sum_i \alpha_i\cdot g(\sigma_i) + \beta,\ \ \ \beta \in B_q(\mathscr C)\f\]
+- A **co-annotation** associates to a cycle its decomposition in the basis of cohomology generators, up to coboundaries: <br>
+Given \f$\gamma\in Z^q(\mathscr C)\f$, \f$g^*(\gamma)\f$ is a co-annotation of \f$\gamma\f$:
+\f\[\text{if }g^*(\gamma) = \sum_i \alpha_i\cdot \sigma_i,\ \ \ \sigma_i\in C^q\f\]
+\f\[\text{then }\ \ \ \gamma = \sum_i \alpha_i\cdot f^*(\sigma_i) + \beta,\ \ \ \beta \in B^q(\mathscr C)\f\]
+
+Annotations and co-annotations provide convenient tests to decide if to cycles (resp. co-cycles) belong to the same homology (resp. cohomology) class:
+- Two cycles \f$\gamma\f$ and \f$\gamma'\f$ belong to the same homology class if they have the same annotation (ie.\ \f$f(\gamma) = f(\gamma')\f$)
+- Two co-cycles \f$\gamma\f$ and \f$\gamma'\f$ belong to the same cohomology class if they have the same co-annotation (ie.\ \f$g^*(\gamma) = g^*(\gamma')\f$)
+
+Figure \cgalFigureRef{HDVFex_hom_cohom} shows an example of HDVF together with corresponding homology and cohomology generators. In all the package, we use the following colors convention:
+- PRIMARY cells are depicted in green
+- SECONDARY cells are depicted in red
+- CRITICAL cells are depicted in blue
+
+Let us point out that HDVFs can be represented in a "discrete Morse theory" style, with arrows (ie.\ a vector field) pairing a primary and a secondary cell among its co-faces (the vector field associated to a HDVF is not unique). But DGVF are a strict subset of HDVFs:
+- each DGVF is a HDVF (labelling sources of arrows as PRIMARY, targets of arrows as SECONDARY, unpaired cells as CRITICAL always produces an HDVF)
+- but HDVF may produce vector fields with cycles
+
+The following figures show an HDVF and an associated vector field:<br>
+ <img src="HDVF_ex1_A3.png" align="center" width=30%/>
+ <img src="HDVF_ex1_A3_DGVF.png" align="center" width=30%/>
+
+ And the following figure exhibits the vector field representation of a HDVF with a cycle (on the right, the cycle is highlighted in purple):<br>
+  <img src="HDVF_cycle.png" align="center" width=20%/>
+ <img src="HDVF_cycle_highlighted.png" align="center" width=20%/><br>
+ Hence this vector field is a valid HDVF (\f$\partial(S)|_P\f$ is invertible) but it is not a DGVF.
+
+\subsection secHDVFcomput Computing and deforming a HDVF
+
+HDVFs can be computed using previous matrix expression; however, Aldo Gonzalez-Lorenzo has shown that several operations on HDVF can be introduced, not only to built, but also transform HDVFs.
+
+- **A operation** (pairing of two critical cells)<br>
+Given a HDVF \f$X=(P,S,C)\f$, and two CRITICAL cells \f$\gamma_1\in C_q\f$, \f$\gamma_2\in C_{q+1}\f$
+\f\[ A(X) = \Big(P\cup\{\gamma_1\}, S\cup\{\gamma_2\}, C\backslash \{\gamma_1, \gamma_2\}\Big)\f\]
+and the operation is valid (ie.\ produces an HDVF) if and only if:
+\f\[ \langle \mathrm d(\gamma_2),\gamma_1\rangle \text{ is invertible}\f\]
+
+    - Any perfect HDVF can be built by iterating the A operation (with complexity at most \f$\mathscr O(n^3))\f$)
+    - This operation can also complete a HDVF into a perfect HDVF
+
+    Figures below illustrate this operation on a cubical chain complex: (left) an HDVF such that \f$(\gamma_1,\gamma_2)\f$ meet the A condition, (right) the HDVF obtained after \f$A(\gamma_1,\gamma_2)\f$.
+
+<img src="HDVF_ex1_A2_A.png" align="center" width=30%/>
+<img src="HDVF_ex1_A3.png" align="center" width=30%/>
+
+- **R operation** (unpairing of two critical cells)<br>
+Given a HDVF \f$X=(P,S,C)\f$, a PRIMARY cell \f$\pi\in P_q\f$ and a SECONDARY cells \f$\sigma\in S_{q+1}\f$
+\f\[ R(X) = \Big(P\backslash\{\pi\}, S\backslash\{\sigma\}, C\cup \{\pi, \sigma\}\Big)\f\]
+and the operation is valid (ie.\ produces an HDVF) if and only if:
+\f\[ \langle h(\pi),\sigma\rangle \text{ is invertible}\f\]
+
+    Operations A and R are inverse one another.
+
+- **M operation** (exchange a PRIMARY and a CRITICAL cell)<br>
+Given a HDVF \f$X=(P,S,C)\f$, a PRIMARY cell \f$\pi\in P_q\f$ and a CRITICAL cell \f$\gamma_1\in C_{q}\f$
+\f\[ M(X) = \Big(P\backslash \{\pi\} \cup\{\gamma\}, S, C\backslash \{\gamma\}\cup \{\pi\}\Big)\f\]
+and the operation is valid (ie.\ produces an HDVF) if and only if:
+\f\[ \langle \mathrm f(\pi),\gamma\rangle \text{ is invertible}\f\]
+
+    The M operation modifies the homology generator associated to \f$\gamma\f$ (while preserving is associated cohomology generator).
+    The following figure illustrates the M operation.
+
+ \cgalFigureBegin{HDVF_M,HDVF_M.png}
+ M operation on a cubical complex \f$K\f$. **Left:** perfect HDVF over \f$K\f$ with cells \f$(\pi, \gamma)\f$ meeting the M operation condition (blue: homology generator associated with \f$\gamma\f$, pink: cohomology generator associated with \f$\gamma\f$); **Right:** perfect HDVF obtained with \f$M(\pi,\gamma)\f$ and updated generators (associated to \f$\pi\f$ which is now critical).
+ \cgalFigureEnd
+
+ - **W operation** (exchange a SECONDARY and a CRITICAL cell)<br>
+Given a HDVF \f$X=(P,S,C)\f$, a SECONDARY cell \f$\sigma\in S_q\f$ and a CRITICAL cell \f$\gamma_1\in C_{q}\f$
+\f\[ W(X) = \Big(P, S\backslash\{\sigma\}\cup \{\gamma\}, C\backslash \{\gamma\}\cup \{\sigma\}\Big)\f\]
+and the operation is valid (ie.\ produces an HDVF) if and only if:
+\f\[ \langle \mathrm g(\gamma),\sigma\rangle \text{ is invertible}\f\]
+
+    The W operation modifies the cohomology generator associated to \f$\gamma\f$ (while preserving is associated homology generator).
+    The following figure illustrates the W operation.
+
+ \cgalFigureBegin{HDVF_W,HDVF_W.png}
+ W operation on a cubical complex \f$K\f$. **Left:** perfect HDVF over \f$K\f$ with cells \f$(\sigma, \gamma)\f$ meeting the W operation condition (blue: homology generator associated with \f$\gamma\f$, pink: cohomology generator associated with \f$\gamma\f$); **Right:** perfect HDVF obtained with \f$W(\pi,\gamma)\f$ and updated generators (associated to \f$\sigma\f$ which is now critical).
+ \cgalFigureEnd
+
+- **MW operation** (exchange a PRIMARY and a SECONDARY cell)<br>
+    Given a HDVF \f$X=(P,S,C)\f$, a PRIMARY cell \f$\pi\in P_q\f$ and a SECONDARY cell \f$\sigma\in S_{q}\f$
+    \f\[ MW(X) = \Big(P\backslash\{\pi\}\cup \{\sigma\}, S\backslash\{\sigma\}\cup \{\pi\}, C\Big)\f\]
+    and the operation is valid (ie.\ produces an HDVF) if and only if:
+    \f\[ \langle \mathrm dh(\pi),\sigma\rangle \text{ and } \langle h\mathrm d(\sigma),\pi\rangle\text{ are invertible}\f\]
+
+    The MW operation modifies both homology and cohomology generators.
+    The following figure illustrates the W operation.
+
+ \cgalFigureBegin{HDVF_MW,HDVF_MW.png}
+ MW operation on a cubical complex \f$K\f$. **Top:** perfect HDVF over \f$K\f$ with cells \f$(\pi_1, \sigma_1)\f$ and \f$(\pi_2, \sigma_2)\f$ meeting the MW operation condition (blue: homology generators, pink: cohomology generators); **Down, left:** perfect HDVF obtained with \f$MW(\pi_1,\sigma_1)\f$ and updated generators (a cohomology generator is modified); **Down, right:** perfect HDVF obtained with \f$MW(\pi_2,\sigma_2)\f$ and updated generators (a homology generator is modified).
+ \cgalFigureEnd
+
+Recent works by Yann-Situ Gazull [YSG1], [YSG2] proved that (when \f$A\f$ is a field):
+- the space of perfect HDVFs is connected under M, W and MW operations (see [YSG2])
+- homology bases computed by HDVFs and other computational approaches (Smith Normal Form, persistent homology, tripartitions...) can be precisely characterised as **explicit basis** (see [YSG1])
+- any perfect HDVF can be adapted in time \f$\mathscr O(kn^2)\f$ to compute a free set of generators of size \f$k\f$ (e.g. a given basis)
+
+These results will soon be implemented in the HDVF package.
+
+\subsection secHDVFpersistent_homology Persistent homology
+
+\subsection secHDVFreduced_Alexander Relative homology and Alexander duality
+
+\section secHDVFHDVFpackage The HDVF package
+
+This section describes the HDVF package implementing the theory described in previous sections.
+
+In order to compute the homology/cohomology of a discrete object with an HDVF, the process is as follows:
+-# Define the ring of coefficients type
+-# Build a chain complex from one of the models of `AbstractChainComplex` (`Geometric_simplicial_chain_complex`, `Abstract_simplicial_chain_complex` or `Cubical_chain_complex`)
+-# Build a HDVF from this chain complex
+-# Compute a perfect HDVF (fast mode or random mode)
+
+ \cgalFigureBegin{HDVFoverview_HDVF,HDVFoverview_HDVF.png}
+ Overview of homology/cohomology computation with the HDVF package.
+ \cgalFigureEnd
+
+\subsection secHDVFCoefficientType Defining the ring of coefficients
+
+The ring of coefficients must be a model of the `IntegralDomainWithoutDivision` concept:
+- the ring \f$\mathbb Z\f$
+\code
+    typedef int CoefficientType ;
+\endcode
+- a field \f$\mathbb Z/p\mathbb Z\f$
+\code
+    #include <CGAL/HDVF/Zp.h>
+    typedef Zp<p> CoefficientType ;
+\endcode
+with \f$p\f$ a given prime number.
+- the field \f$\mathbb Z/2\mathbb Z\f$ (the class`Zp` can also be used, but `Z2`is optimized for \f$p=2\f$)
+\code
+    #include <CGAL/HDVF/Z2.h>
+    typedef Z2 CoefficientType ;
+\endcode
+
+\subsection secHDVFChainComplex Building a chain complex
+
+The HDVF package provides several chain complexes: `Abstract_simplicial_chain_complex`, `Geometric_simplicial_chain_complex` and `Cubical_chain_complex`.
+
+In order to build such a chain complex, discrete objects are first imported into a temporary combinatorial structure storing input cells. The chain complex is then built from this intermediate structure by adding all the faces of cells and computing the corresponding boundary matrix.
+
+Let us point out that whatever the type of chain complex, cells are indexed along each dimension; hence a cell is identified by an index and a dimension.
+
+- `Abstract_simplicial_chain_complex` (see section \ref secHDVFsimplicial_chain_complex) are built using the `Mesh_object_io` intermediate class from abstract simplices (vertices have no coordinates). Therefore, such complexes are usually computed from `.simp` files (a simple file format enumerating on each line vertices of simplices, see section secHDVFsimpFile):
+    \code
+        using Complex = Abstract_simplicial_chain_complex<CoefficientType> ;
+        Mesh_object_io mesh ;
+        mesh.read_simp(filename);
+        Complex complex(mesh);
+    \endcode
+- `Simplicial_chain_complex` are also built using the `Mesh_object_io` intermediate structure to import mesh files (`.off` format):
+    \code
+        using Complex = Simplicial_chain_complex<CoefficientType> ;
+        Mesh_object_io mesh ;
+        mesh.read_off(filename);
+        Complex complex(mesh);
+    \endcode
+- `Cubical_chain_complex` are built using the `Cub_object_io` intermediate structure by importing either simple `.cub` files (enumerating Khalimsky coordinates of cells, see section \ref secHDVFcubFile) or `.pgm` files.
+    - For `.cub` files enumerating Khalimsky coordinates, cells are imported and their faces are added to produce a proper complex.
+        \code
+            using Complex = Cubical_chain_complex<CoefficientType> ;
+            Cub_object_io mesh ;
+            mesh.read_cub(filename, true);
+            Complex complex(mesh);
+        \endcode
+    - For `.pgm` files (ie.\ for binary images or volumes), two constructions are available for cubical complexes:
+        - PRIMAL, encoding \f$3^d-1\f$-adjacency  (where cells are imported as such, see figure \cgalFigureRef{HDVF_primal_dual}, middle)
+        \code
+            using Complex = Cubical_chain_complex<CoefficientType> ;
+            Cub_object_io mesh ;
+            // Read pgm file (store Khalimsky coordinates)
+            mesh.read_pgm(filename, true);
+            // Build the PRIMAL complex associated
+            Complex complex(mesh, CGAL::Homological_discrete_vector_field::Complex::PRIMAL);
+        \endcode
+        - DUAL, encoding \f$2d\f$-adjacency (where cubes produce vertices... see figure \cgalFigureRef{HDVF_primal_dual} right)
+        \code
+            using Complex = Cubical_chain_complex<CoefficientType> ;
+            Cub_object_io mesh ;
+            // Read pgm file (store voxel coordinates)
+            mesh.read_pgm(filename, false);
+            // Build the DUAL complex associated
+            Complex complex(mesh, CGAL::Homological_discrete_vector_field::Complex::DUAL);
+        \endcode
+
+\cgalFigureBegin{HDVF_primal_dual,HDVF_primal_dual.png}
+Primal and dual constructions of cubical chain complexes from a binary volume. **Left:** a 3D binary volume; **Middle:** PRIMAL cubical (chain) complex associated to the volume; **Right:** DUAL cubical (chain) complex associated to the volume).<br>
+    Figure taken from [AGL, 2016].
+\cgalFigureEnd
+
+\subsubsection secHDVFsimpFile The `.simp` file format
+
+The `simp` file format is a very simple text format listing, encoding simplicial complexes. It enumerates simplices line by line (described by the indices of their vertices).
+The file doesn't need to provide an exhaustive enumeration of simplices (with all their faces), indeed, the simplicial chain complex constructors fill them with faces. Figure \cgalFigureRef{HDVF_simp_file} illustrate the format on a simple simplicial complex.
+
+\cgalFigureBegin{HDVF_simp_file,HDVF_simp_file.png}
+Example of `simp` file. **Left:** a simplicial complex; **Right:** associated `simp` files describing maximal simplices.
+\cgalFigureEnd
+
+\subsubsection secHDVFcubFile The `.cub` file format
+
+The `cub` file format is a very simple text format encoding binary volumes:
+- first line: dimension
+- second line: size of Khalimsky coordinates of the bounding box along each axis
+- following lines: list of cells in Khalimsky coordinates
+
+The file doesn't need to provide an exhaustive enumeration of cells (with all their faces), indeed, the cubical chain complex constructors fill them with faces. Figure \cgalFigureRef{HDVF_cub_file} illustrate the format on a simple cubical complex.
+
+\cgalFigureBegin{HDVF_cub_file,HDVF_cub_file.png}
+Example of `cub` file. **Left:** a cubical complex; **Right:** associated `cub` files describing maximal cells.
+\cgalFigureEnd
+
+\subsection secHDVFbuild Building and computing a HDVF
+
+An empty HDVF is built from a chain complex:
+\code
+    using HDVF_type = Hdvf<CoefficientType, Complex> ;
+    HDVF_type hdvf(complex, OPT_FULL) ;
+\endcode
+
+Several HDVF options are available:
+- `OPT_D`: fastest option - compute only the reduced boundary \f$\mathrm d\f$ (under this option, HDVFs provide only low level homological information: Betti numbers and critical cells
+- `OPT_F`, `OPT_G`: compute the reduction partially (only \f$f\f$ for cohomology generators or \f$g\f$ for homology generators)
+- `OPT_FULL`: compute the full reduction
+
+\subsubsection secHDVFcomputePerfect Computing a perfect HDVF automatically
+
+The HDVF obtained is empty (all cells are critical). In order to compute homology/cohomology that is a perfect HDVF, the package provides two methods:
+- `compute_perfect_hdvf` (fastest method:  cells are paired using the A operation by decreasing dimension using the first valid pair)
+- `compute_rand_perfect_hdvf` (slowest: compute all valid pairs for the A operation are select one randomly)
+
+\code
+    hdvf.compute_perfect_hdvf();
+\endcode
+or
+\code
+    hdvf.compute_rand_perfect_hdvf();
+\endcode
+
+When the ring of coefficients is a field, these methods always produce a perfect HDVF. However for \f$\mathbb Z\f$-(co)homology (without torsion), the existence of complexes admitting no perfect HDVF is still an open question (and the answer is negative for "reasonable" objects, indeed we still have no counter example).
+The method `is_perfect_hdvf` checks if a HDVF is perfect.
+
+\subsubsection secHDVFcomputePerfectA Computing a perfect HDVF with A operations
+
+For specific applications, HDVFs can be built by iterating the A operation:
+\code
+    hdvf.A(id1, id2, q) ;
+\endcode
+performs a pairing between the cells \f$\gamma_1\f$ (index `id1` in dimension `q`) and \f$\gamma_2\f$ (index `id2` in dimension `q+1`).
+
+The package provide a family of methods to search for valid pairings:
+\code
+    Cell_pair pair;
+    std::vector<Cell_pair> pairs;
+    bool found;
+
+    // Search for the first valid pair in dimension `q`
+    pair = find_pair_A(q, found);
+    // Search for the first cell (of dimension `q-1` or `q+1`) forming a valid pair with tau (index `id`, dimension `q`)
+    pair = find_pair_A(q, found, id);
+    // Search for all valid pairs in dimension `q`
+    pairs = find_pairs_A(q, found);
+    // Search for all cells (of dimension `q-1` or `q+1`) forming a valid pair with tau (index `id`, dimension `q`)
+    pairs = find_pairs_A(q, found, id);
+\endcode
+
+\subsection secHDVFoutput Extracting and saving homological/cohomological information
+
+\warning All `vtk` exports are only available for geometric chain complexes (ie.\ `Simplicial_chain_complex` and `Cubical_chain_complex`).
+
+How to extract (co)homological information then? Various solutions and output are available:
+- **HDVF extraction / visualization**:
+    The method `psc_labels` extracts cells flags in any dimensions. The complex annotated with labels can then be exported to `vtk` for visualization. The `vtk` output includes a `Label` property encoding the HDVF PSC_flag (-1 for PRIMARY, 1 for SECONDARY and 0 for CRITICAL cells), and a `CellId` property storing cells indices:
+    \code
+        // Labels extraction
+        vector<vector<int> > labels = hdvf.psc_labels() ;
+        // VTK export
+        Complex::chain_complex_to_vtk(complex, filename, &labels) ;
+    \endcode
+- **Betti numbers**:
+    For a perfect HDVF, the Betti number in dimension \f$q\f$ equals \f$|C_q|\f$, the number of critical cells.
+    Critical cells can be obtained:
+    - either in all dimensions:
+    \code
+        std::vector<std::vector<size_t> > criticals =  hdvf.psc_flags(CGAL::Homological_discrete_vector_field::CRITICAL);
+    \endcode
+    - or in a given dimension `q`:
+    \code
+        std::vector<size_t> criticals_q = hdvf.psc_flags(CGAL::Homological_discrete_vector_field::CRITICAL, q) ;
+    \endcode
+- **Homology/cohomology generators**
+    Each critical cell \f$\gamma\f$ bears an associated homology generator: \f$g(\gamma)\f$ and cohomology generator: \f$f^*(\gamma)\f$.
+    - Generators associated to \f$\gamma\f$ (index `id`, dimension `q`) can be extracted as column chains (ie.\ column `Sparse_chain`):
+    \code
+        HDVF_type::Column_chain hom_generator(hdvf.homology_chain(id, q)) ;
+        HDVF_type::Column_chain cohom_generator(hdvf.cohomology_chain(id, q)) ;
+    \endcode
+    These chains can be visited as described in section \ref secHDVFmatrixChain.
+    - The HDVF package also provides a vtk exporter for such chains (cells with a non zero coefficient are exported). For instance:
+    \code
+        Complex::chain_to_vtk(complex, filename, hom_generator, q, id) ;
+    \endcode
+    If `id` is provided, a `Label` integer property is added to the vtk file (0 for the cell `id`, 2 for other cells).
+- **Reduction maps**
+    Matrices of the reduction (\f$F\f$, \f$G\f$, \f$H\f$ and \f$D\f$) can be also obtained using `matrix_f`, `matrix_g`, `matrix_h` and `matrix_dd`.
+- **Automated vtk export**
+    In order to ease `vtk` export, an "automated" method `CGAL::IO::write_VTK` is also provided, exporting in distinct files, the HDVF flags and all the generators computed (homology / cohomology depending of the HDVF option):
+    \code
+            CGAL::IO::write_VTK<CoefficientType, Complex>(hdvf, complex, outfile_root, co_faces) ;
+    \endcode
+    where `outfile_root` is the root of the `.vtk` file names and `co_faces` is a %Boolean (if `true`, cohomology generators are represented through their cofaces).
+- HDVFs and reductions can also be saved / loaded to `.hdvf` files (a simple text file format described in section \ref secHDVFHDVF_format) with `insert_hdvf_reduction()` and `extract_hdvf_reduction()`.
+
+\subsubsection secHDVFHDVF_format The `.hdvf` file format
+
+The `.hdvf` file format is a simple text format used to save HDVFs. It relies on the `.osm` format described in section \ref secHDVFOSM_format.
+
+\warning The file format does not store the underlying complex, therefore the dimension and the number of cells in each dimension are used to check the coherence of the `.hdvf` file with the underlying complex provided by the user.
+
+- HDVF type (0 for HDVF and reduction, 1 for HDVF only)
+- Dimension of the complex
+- List of the number of cells in each dimension (one dimension by row)
+- HDVF option (OPT_BND | OPT_F | OPT_G | OPT_FULL)
+- List of the flags of cells in each dimension (one dimension by row: -1 for PRIMARY, 1 for SECONDARY, 0 for CRITICAL)
+- According to HDVF option: output F matrices in each dimension (`.osm` format)
+- According to HDVF option: output G matrices in each dimension (`.osm` format)
+- Output H matrices in each dimension (`.osm` format)
+- Output DD matrices in each dimension (`.osm` format)
+
+\subsection secHDVFexHDVF HDVF examples
+
+The following examples exhibit the computation of a perfect HDVF over a simplicial (resp. cubical) chain complex built from a `.off` file (resp. a `.pgm` file). Generators are then exported to `.vtk`.
+\cgalExample{HDVF/hdvf_simplicial.cpp}
+\cgalExample{HDVF/hdvf_cubical.cpp}
+
+\subsection secHDVFmatrixChain A few words about sparse matrices and chains
+
+All maps (boundary operator, \f$f\f$, \f$g\f$, \f$h\f$ and \f$\mathrm d\f$) are stored in sparse matrices and chains are stored in sparse chains.
+The `OSM` library provides an implementation of the `SparseMatrix` and `SparseChain` concepts, that is sparse structures optimized for topological operations (block operations).
+
+The structure is based on a mapped matrix implementation, hence, matrices and chains are either column or row major (`CGAL::OSM::COLUMN` or `CGAL::OSM::ROW` type).
+
+In order to speed up both block operations and iterators (over non empty columns or row according to the type), a **bitboard** structure is used (inherited from chess programming community). Hence, in order to iterate over the major dimension of a matrix, one should actually iterate a bitboard (providing indices of non empty columns / rows). See `Sparse_matrix` documentation for more details.
+
+Iterating over sparse chains them comes to iterate over a map.
+
+The following example builds a column major matrix, iterate over its non empty columns (ie.\ chains) and then iterate over each chain.
+
+\cgalExample{HDVF/matrix_chain.cpp}
+
+All standard linear algebra operators are provided (and optimized according to the matrix/chain types).
+
+\subsubsection secHDVFOSM_format The `.osm` file format
+
+The `.osm` file format is a simple text format to store OSM sparse matrices:
+- Matrix type (0 for COLUMN, 1 for ROW)
+- Number of rows and number of columns (separated by a space)
+- Number of coefficients
+- List of coefficients (one by row):<br>
+    `id_row  id_col  value`
+
+\subsection secHDVFoperations HDVF operations
+
+\subsection secHDVFannotations HDVFs and annotations
+
+\subsection secHDVFpersistentHom_package HDVFs and persistent homology
+
+\subsection secHDVFAlexander_package HDVFs, Alexander duality and relative homology
+
+
+
+ [AGL, 2017] Aldo Gonzalez-Lorenzo, Alexandra Bac, Jean-Luc Mari, Pedro Real. Allowing cycles in discrete Morse theory, Topology and its Applications, Volume 228, 2017, Pages 1-35.
+
+ [AGL, 2016] Aldo Gonzalez-Lorenzo. Computational Homology Applied to Discrete Objects. Discrete Mathematics [PhD]. Aix-Marseille Université; Universidad de Sevilla, 2016.
+
+ [AGL, 2025] Aldo Gonzalez-Lorenzo, Alexandra Bac and Yann-Situ Gazull. A constructive approach of Alexander duality. J Appl. and Comput. Topology 9, 2 (2025).
+
+ [YSG1,  2025] Yann-Situ Gazull, Aldo Gonzalez-Lorenzo and Alexandra Bac. Characterization of the computed homology and cohomology bases. DGMM 2025 (to appear).
+ [YSG1,  2025] Yann-Situ Gazull, Aldo Gonzalez-Lorenzo and Alexandra Bac. Characterization of the computed homology and cohomology bases. DGMM 2025 (to appear).
+
+ [YSG2,  2025] Yann-Situ Gazull, Aldo-Gonzalez-Lorenzo and Alexandra Bac. Space of homological computations: connectivity and paths (submitted to Journal of Applied and Computational Topology).
+
+ [Munkres] J. Munkres. Elements of Algebraic Topology, Addison Wesley Publishing Company, (1984)
+
+*/
+} /* namespace CGAL */

HDVF/doc/HDVF/PackageDescription.txt

@@ -0,0 +1,47 @@
+// PRETTY HDVF should equal the project title in Doxyfile.in
+
+/// \defgroup PkgHDVFRef Homological Discrete Vector Fields Reference
+/// \defgroup PkgHDVFConcepts Concepts
+/// \ingroup PkgHDVFRef
+
+/// \defgroup PkgHDVFAlgorithmClasses Sparse Matrices and Vectors
+/// \ingroup PkgHDVFRef
+
+/// \defgroup PkgHDVFTraitsClasses Traits Classes
+/// \ingroup PkgHDVFRef
+
+/// \defgroup PkgHDVFIOClasses IO Classes
+/// \ingroup PkgHDVFRef
+
+
+/*!
+\addtogroup PkgHDVFRef
+\cgalPkgDescriptionBegin{Homological Discrete Vector Fields, PkgHDVF}
+\cgalPkgPicture{HDVF-small.png}
+
+\cgalPkgSummaryBegin
+\cgalPkgAuthors{Alexandra Bac}
+\cgalPkgDesc{
+The package provides an implementation of Homological Discrete %Vector Fields, a framework computing homology and cohomology information (Betti numbers, generators, annotation functions...) for discrete objects (including persistent homology, Alexander duality isomorphism and reduced homology).}
+\cgalPkgManuals{Chapter_HDVF,PkgHDVFRef}
+\cgalPkgSummaryEnd
+
+\cgalPkgShortInfoBegin
+\cgalPkgSince{6.2}
+\cgalPkgBib{cgal:b-hdvf}
+\cgalPkgLicense{\ref licensesGPL "GPL"}
+
+\cgalPkgShortInfoEnd
+
+\cgalPkgDescriptionEnd
+
+\cgalClassifedRefPages
+
+\cgalCRPSection{Concepts}
+- `AbstractChainComplex`
+
+\cgalCRPSection{Classes}
+- `CGAL::Homological_discrete_vector_field::Abstract_simplicial_chain_complex`
+
+
+*/

HDVF/doc/HDVF/dependencies

@@ -0,0 +1,9 @@
+Manual
+Kernel_23
+Kernel_d
+STL_Extension
+Algebraic_foundations
+Number_types
+Circulator
+Stream_support
+Surface_mesh

HDVF/doc/HDVF/examples.txt

@@ -0,0 +1,9 @@
+/*!
+\example HDVF/main_dual_hdvf.cpp
+\example HDVF/main_hdvf.cpp
+\example HDVF/main_per_hdvf.cpp
+\example HDVF/hdvf_simplicial.cpp
+\example HDVF/hdvf_surface_mesh_simplicial.cpp
+\example HDVF/hdvf_cubical.cpp
+\example HDVF/matrix_chain.cpp
+*/

HDVF/examples/HDVF/CMakeLists.txt

@@ -0,0 +1,26 @@
+# Created by the script cgal_create_cmake_script
+# This is the CMake script for compiling a CGAL application.
+
+cmake_minimum_required(VERSION 3.12...3.31)
+project(HDVF_Examples)
+
+find_package(CGAL REQUIRED)
+
+
+find_package(Eigen3 3.1.0 QUIET) #(requires 3.1.0 or greater)
+include(CGAL_Eigen3_support)
+
+
+
+# create a target per cppfile
+file(
+  GLOB cppfiles
+  RELATIVE ${CMAKE_CURRENT_SOURCE_DIR}
+  ${CMAKE_CURRENT_SOURCE_DIR}/*.cpp)
+foreach(cppfile ${cppfiles})
+  create_single_source_cgal_program("${cppfile}")
+endforeach()
+
+if(TARGET CGAL::Eigen3_support)
+target_link_libraries(hdvf_simplicial_4 PUBLIC CGAL::Eigen3_support)
+endif()

HDVF/examples/HDVF/arguments.h

@@ -0,0 +1,394 @@
+#ifndef ARGUMENTS_H
+#define ARGUMENTS_H
+
+#include <CGAL/HDVF/Hdvf.h>
+
+using namespace std;
+using namespace CGAL;
+using namespace Homological_discrete_vector_field;
+
+enum class Algorithm { HDVF, DualHDVF, PerHDVF };
+enum class InputFormat { SIMP, OFF, CUB, PGM };
+enum class StarFiltrStd { FiltrX, FiltrY, FiltrZ };
+
+struct Options
+{
+    std::string outfile_root = "res", in_file ;
+    Algorithm algorithm = Algorithm::HDVF ;
+    InputFormat in_format = InputFormat::OFF ;
+    int HDVF_opt = CGAL::Homological_discrete_vector_field::OPT_FULL ;
+    int scalar = 0 ; // 0 : Z, n : Z/nZ
+    bool with_export = true ; // Export reduction information to a file
+    bool with_vtk_export = false ; // Export generators/cogenerators and PSC labels as vtk
+    bool with_output = true ; // Consol output
+    bool with_frame = true ; // Only for persistent homology (cub and pgm in_files)
+    bool loop = false ; // Only for HDVF (enter an interaction loop for M, W, MW)
+    double BB_ratio = 1.5 ; // Only for persistent homolgy (off in_files: size of the bounding sphere)
+    bool primal = true ; // Only for cub_complexes (for dual construction, set primal to false)
+    StarFiltrStd star_filtr = StarFiltrStd::FiltrZ ; // Only for per_hdvf (when filtrataion is a standard lower star filtration)
+    bool random = false ; // Activate the computation of random perfect HDVF
+    bool verbose = false ; // Activate verbose mode (matrix output after each A operation for HDVF computation)
+    bool co_faces = false ; // Output co-faces of cohomology generators (or not)
+
+    std::string comment() const
+    {
+        std::string str;
+        if (algorithm == Algorithm::HDVF)
+            str = "HDVF computation";
+        else if (algorithm == Algorithm::DualHDVF)
+            str = "Alexander duality HDVF computation";
+        else if (algorithm == Algorithm::PerHDVF)
+            str = "Persistent HDVF computation";
+        else
+            str = "undefined";
+        return str;
+    }
+};
+
+inline void usage() {
+    cout << "usage:" << endl ;
+    cout << "\talgorithm [options] in_file" << endl ;
+    cout << "with:" << endl ;
+    cout << "algorithm = hdvf (for hdvf (co)homology compuation)" << endl ;
+    cout << "\t dual_hdvf (for hdvf Alexander duality)" << endl ;
+    cout << "\t per_hdvf (for hdvf persistent (co)homology)" << endl ;
+    cout << "and options:" << endl ;
+    cout << "\t-outfile_root = xxx : sets the outfile root for outfiles names to xxx (default: res)" << endl ;
+    cout << "\t-opt (bnd|f|g|full) : HDVF option (default: full)" << endl ;
+    cout << "\t-s n : ring used to compute homology, 0: Z, n: Z/nZ (default: 0)" << endl ;
+    cout << "\t-rand : compute a random perfect HDVF (default: false, computing a random HDVF is slower)" << endl ;
+    cout << "\t-verbose : verbose HDVF computation, output reduction matrices after each A (default: false)" << endl ;
+    cout << "\t-with_output / -no_output : output / do not output the HDVF and reduction to the console (default: with_output)" << endl ;
+    cout << "\t-with_export / -no_export : export / do not export the HDVF and reduction to a file (default: with_export)" << endl ;
+    cout << "\t-with_vtk_export / -no_vtk_export : export / do not export the HDVF and reduction to vtk (default: no_export)" << endl ;
+    cout << "\t-cofaces_cohomology / -no_cofaces_cohomology : do not export co-faces of cohomology generators (default: co_faces_cohomology)" << endl ;
+    cout << "For CubComplex only:" << endl ;
+    cout << "\t-primal / -dual: build the primal / dual complex associated to data (default: primal)" << endl ;
+    cout << "For Alexander duality only:" << endl ;
+    cout << "\tFor SimpComplex:" << endl ;
+    cout << "\t\t-BB_ratio x : ratio for the closing bounding sphere (default: 1.5)" << endl ;
+    cout << "\tFor CubComplex:" << endl ;
+    cout << "\t\t-with_frame / -no_frame : enclose / do not inclose the complex in a larger bounding box (default: with_frame, adds 1 pixel around)" << endl ;
+    cout << "For persistent HDVF only:" << endl ;
+    cout << "\t-lower (x|y|z) : compute lower star filtration with degree x, y or z" << endl ;
+    cout << "\tFor other filtrations, define the filtration into the main and set #define OwnFiltration" << endl;
+}
+
+inline Options read_arguments_hdvf(int argc, char** argv) {
+    Options opt_res ;
+
+    // Minumum number of arguments: 2
+    // HDVF_algorithm file
+
+    if (argc <= 2)
+    {
+        cout << "HDVF wrong number of arguments. HDVF -h for help" << endl ;
+        throw "HDVF wrong number of arguments" ;
+    }
+    // Read algorithm from argv[0] (last word of the path - separated by /)
+    const string algo(argv[0]) ;
+    string word1 ;
+    {
+        std::stringstream string_parse(algo);
+        std::string tmp_word;
+
+        std::getline(string_parse, tmp_word, '/') ;
+        word1 = tmp_word ;
+        while(std::getline(string_parse, tmp_word, '/'))
+        {
+            word1 = tmp_word ;
+        }
+    }
+
+    if (word1.compare("main_hdvf") == 0)
+        opt_res.algorithm = Algorithm::HDVF ;
+    else if (word1.compare("main_dual_hdvf") == 0)
+        opt_res.algorithm = Algorithm::DualHDVF ;
+    else if (word1.compare("main_per_hdvf") == 0)
+        opt_res.algorithm = Algorithm::PerHDVF ;
+    else
+    {
+        cerr << "HDVF algorithm unknown. HDVF -h for help" << endl ;
+        throw "HDVF algorithm unknown" ;
+    }
+
+    // Read in_file from argv[argc-1] ;
+    opt_res.in_file = argv[argc-1] ;
+    // Deduce the file type from the file name
+    string word ;
+    {
+        // Get the file name (last word of the path - separated by /)
+        std::stringstream string_parse(opt_res.in_file);
+        std::string tmp_word;
+
+        std::getline(string_parse, tmp_word, '/') ;
+        word = tmp_word ;
+        while(std::getline(string_parse, tmp_word, '/'))
+        {
+            word = tmp_word ;
+        }
+    }
+    {
+        // Get the extension (last word of the word - separated by .)
+        string filename = word ;
+        std::stringstream string_parse(filename);
+        std::string tmp_word;
+
+        std::getline(string_parse, tmp_word, '.') ;
+        word = tmp_word ;
+        while(std::getline(string_parse, tmp_word, '.'))
+        {
+            word = tmp_word ;
+        }
+    }
+    // Assign in_format accordingly
+    if (word.compare("simp")==0)
+        opt_res.in_format = InputFormat::SIMP ;
+    else if (word.compare("off")==0)
+        opt_res.in_format = InputFormat::OFF ;
+    else if (word.compare("cub")==0)
+        opt_res.in_format = InputFormat::CUB ;
+    else if (word.compare("pgm")==0)
+        opt_res.in_format = InputFormat::PGM ;
+    else
+    {
+        cerr << "HDVF input format unkown. HDVF -h for help" << endl ;
+        throw "HDVF input format unknown" ;
+    }
+
+
+    // Read next options
+    for(int i = 1; i<argc-1; ++i)
+    {
+        string tmp(argv[i]) ;
+        if (tmp.compare("-h")==0)
+        {
+            usage() ;
+        }
+        else if (tmp.compare("-opt")==0) // HDVF option
+        {
+            if (i == argc-2)
+            {
+                cerr << "HDVF: -opt option requires an option. HDVF -h for help" << endl ;
+                throw "HDVF: -opt option requires an option" ;
+            }
+            string tmp2(argv[++i]) ;
+            if (tmp2.compare("bnd")==0)
+                opt_res.HDVF_opt = OPT_BND ;
+            else if (tmp2.compare("g")==0)
+                opt_res.HDVF_opt = OPT_G ;
+            else if (tmp2.compare("f")==0)
+                opt_res.HDVF_opt = OPT_F ;
+            else if (tmp2.compare("full")==0)
+                opt_res.HDVF_opt = OPT_FULL ;
+            else
+            {
+                cerr << "HDVF option unkown. HDVF -h for help" << endl ;
+                throw "HDVF option unknown" ;
+            }
+        }
+        else if (tmp.compare("-s")==0) // scalar
+        {
+            if (i == argc-2)
+            {
+                cerr << "HDVF: -s option requires an integer. HDVF -h for help" << endl ;
+                throw "HDVF: -s option requires an integer" ;
+            }
+            string tmp2(argv[++i]) ;
+            opt_res.scalar = stoi(tmp2) ;
+        }
+        else if (tmp.compare("-rand")==0) // random HDVF
+        {
+            opt_res.random = true ;
+        }
+        else if (tmp.compare("-verbose")==0) // verbose HDVF computation
+        {
+            opt_res.verbose = true ;
+        }
+        else if (tmp.compare("-BB_ratio")==0) // scalar
+        {
+            if (i == argc-2)
+            {
+                cerr << "HDVF: -BB_ratio option requires a float. HDVF -h for help" << endl ;
+                throw "HDVF: -BB_ratio option requires a float" ;
+            }
+            else if ((opt_res.algorithm != Algorithm::PerHDVF) || (opt_res.in_format != InputFormat::OFF))
+            {
+                cerr << "HDVF: -BB_ratio only for per_HDVF with off input. HDVF -h for help" << endl ;
+                throw "HDVF: -BB_ratio only for per_HDVF with off input" ;
+            }
+            string tmp2(argv[++i]) ;
+            opt_res.BB_ratio = stod(tmp2) ;
+        }
+        else if (tmp.compare("-with_frame")==0) // scalar
+        {
+            if (! ((opt_res.in_format == InputFormat::CUB) || (opt_res.in_format == InputFormat::PGM)))
+            {
+                cerr << "HDVF: -with_frame only for cub or pgm input. HDVF -h for help" << endl ;
+                throw "HDVF: -with_frame only for cub or pgm input" ;
+            }
+            opt_res.with_frame = true ;
+        }
+        else if (tmp.compare("-no_frame")==0) // scalar
+        {
+            if (! ((opt_res.in_format == InputFormat::CUB) || (opt_res.in_format == InputFormat::PGM)))
+            {
+                cerr << "HDVF: -no_frame only for cub or pgm input. HDVF -h for help" << endl ;
+                throw "HDVF: -no_frame only for cub or pgm input" ;
+            }
+            opt_res.with_frame = false ;
+        }
+        else if (tmp.compare("-with_output")==0) // scalar
+        {
+            opt_res.with_output = true ;
+        }
+        else if (tmp.compare("-no_output")==0) // scalar
+        {
+            opt_res.with_output = false ;
+        }
+        else if (tmp.compare("-with_export")==0) // scalar
+        {
+            opt_res.with_export = true ;
+        }
+        else if (tmp.compare("-no_export")==0) // scalar
+        {
+            opt_res.with_export = false ;
+        }
+        else if (tmp.compare("-with_vtk_export")==0) // scalar
+        {
+            if (opt_res.in_format == InputFormat::SIMP)
+            {
+                cerr << "HDVF: -with_vtk_export requires a geometry. HDVF -h for help" << endl ;
+                throw "HDVF: -with_vtk_export requires a geometry" ;
+            }
+            opt_res.with_vtk_export = true ;
+        }
+        else if (tmp.compare("-no_vtk_export")==0) // scalar
+        {
+            opt_res.with_vtk_export = false ;
+        }
+        else if (tmp.compare("-cofaces_cohomology")==0) // scalar
+        {
+            opt_res.co_faces = true ;
+        }
+        else if (tmp.compare("-no_cofaces_cohomology")==0) // scalar
+        {
+            opt_res.co_faces = false ;
+        }
+        else if (tmp.compare("-loop")==0) // scalar
+        {
+            if (!((opt_res.algorithm == Algorithm::HDVF) && (opt_res.HDVF_opt & OPT_FULL)))
+            {
+                cerr << "HDVF: -loop available only for HDVF computation and OPT_FULL required. HDVF -h for help" << endl ;
+                throw "HDVF: -loop available only for HDVF computation and OPT_FULL required" ;
+            }
+            opt_res.loop = true ;
+        }
+        else if (tmp.compare("-primal")==0) // scalar
+        {
+            if ((opt_res.in_format != InputFormat::CUB) && (opt_res.in_format != InputFormat::PGM))
+            {
+                cerr << "HDVF: -primal available only for cubical data. HDVF -h for help" << endl ;
+                throw "HDVF: -primal available only for cubical data" ;
+            }
+            opt_res.primal = true ;
+        }
+        else if (tmp.compare("-dual")==0) // scalar
+        {
+            if ((opt_res.in_format != InputFormat::CUB) && (opt_res.in_format != InputFormat::PGM))
+            {
+                cerr << "HDVF: -dual available only for cubical data. HDVF -h for help" << endl ;
+                throw "HDVF: -dual available only for cubical data" ;
+            }
+            opt_res.primal = false ;
+        }
+        else if (tmp.compare("-lower")==0) // scalar
+        {
+            if (opt_res.algorithm != Algorithm::PerHDVF)
+            {
+                cerr << "HDVF: -lower available only for perHDVF for lower star filtration. HDVF -h for help" << endl ;
+                throw "HDVF: -lower available only for perHDVF for lower star filtration" ;
+            }
+            if (i == argc-2)
+            {
+                cerr << "HDVF: -lower needs an argument x, y or z. HDVF -h for help" << endl ;
+                throw "HDVF: -lower needs an argument x, y or z" ;
+            }
+            else
+            {
+                string tmp(argv[++i]) ;
+
+                if (tmp.compare("x")==0)
+                    opt_res.star_filtr = StarFiltrStd::FiltrX ;
+                else if (tmp.compare("y")==0)
+                    opt_res.star_filtr = StarFiltrStd::FiltrY ;
+                else if (tmp.compare("z")==0)
+                    opt_res.star_filtr = StarFiltrStd::FiltrZ ;
+                else
+                {
+                    cerr << "HDVF: -lower needs an argument x, y or z. HDVF -h for help" << endl ;
+                    throw "HDVF: -lower needs an argument x, y or z" ;
+                }
+            }
+        }
+        else // unknown option
+        {
+            string tmp("HDVF: ") ;
+            tmp += argv[i] ;
+            tmp += " unknown option. HDVF -f for help" ;
+            cerr << tmp << ". HDVF -h for help" << endl ;
+            throw tmp ;
+        }
+    }
+    return opt_res ;
+}
+inline std::ostream& operator<< (std::ostream& out, const Options& opt)
+{
+    out << "outfile_root: " << opt.outfile_root << endl ;
+
+    out << "in_file: " << opt.in_file << endl ;
+
+    out << "algorithm: " << opt.comment() << endl ;
+
+    out << "in_format: " ;
+    if (opt.in_format == InputFormat::OFF)
+        out << "OFF" ;
+    else if (opt.in_format == InputFormat::SIMP)
+        out << "SIMP" ;
+    else if (opt.in_format == InputFormat::CUB)
+        out << "CUB" ;
+    else if (opt.in_format == InputFormat::PGM)
+        out << "PGM" ;
+    out << endl ;
+
+    out << "HDVF_opt: " ;
+    if (opt.HDVF_opt == OPT_FULL)
+        out << "OPT_FULL" ;
+    else if (opt.HDVF_opt == OPT_G)
+        out << "OPT_G" ;
+    else if (opt.HDVF_opt == OPT_F)
+        out << "OPT_F" ;
+    else if (opt.HDVF_opt == OPT_BND)
+        out << "OPT_BND" ;
+    out << endl ;
+
+    out << "scalar: " << opt.scalar << endl ;
+    out << "with_export: " << opt.with_export << endl ; // Export reduction information to a file
+    out << "with_vtk_export: " << opt.with_vtk_export << endl ; // Export generators/cogenerators and PSC labels as vtk
+    out << "co_faces: " << opt.co_faces << endl ; // Export cofaces of cogenerators
+    out << "with_output: " << opt.with_output << endl ; // Consol output
+    out << "with_frame: " << opt.with_frame << endl ; // Only for persistent homology (cub and pgm in_files)
+    out << "BB_ratio: " << opt.BB_ratio << endl ; // Only for persistent homolgy (off in_files: size of the bounding sphere)
+    out << "primal: " << opt.primal << endl ; // Only for cubical data (primal / dual = false)
+    out << "start_filtr: " ;
+    if (opt.star_filtr == StarFiltrStd::FiltrX)
+        out << "x" << endl ;
+    else if (opt.star_filtr == StarFiltrStd::FiltrY)
+        out << "y" << endl ;
+    else if (opt.star_filtr == StarFiltrStd::FiltrZ)
+        out << "z" << endl ;
+    return out ;
+}
+
+#endif // ARGUMENTS_H
+

HDVF/examples/HDVF/data/data_cubical/Amphora_20.pgm

@@ -0,0 +1,159 @@
+P2
+22 13 12
+255
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 1 0 0 0 0 0 1 1 0 0 0 0 0 0 0 
+0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 
+0 0 0 0 0 0 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 1 0 0 0 0 0 0 
+0 0 0 0 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 1 1 1 0 0 1 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 1 1 1 0 1 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 1 1 1 1 1 1 1 0 0 0 0 0 1 0 0 0 0 0 0 
+0 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 
+0 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 
+0 0 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 
+0 0 0 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 1 1 1 0 0 1 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 1 0 0 
+0 0 1 1 1 1 1 1 1 1 1 0 0 0 1 1 1 1 0 1 0 0 
+0 1 1 1 1 1 1 1 1 1 1 1 0 0 0 1 0 0 0 0 0 0 
+0 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 
+0 0 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 
+0 0 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 
+0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 1 0 0 0 0 0 0 
+0 0 0 0 0 1 1 1 1 1 0 0 0 0 0 1 0 0 0 0 0 0 
+0 0 0 1 1 1 1 1 1 1 1 0 0 0 1 0 0 0 0 0 1 0 
+0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 1 0 
+0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 
+0 1 1 1 1 1 1 1 1 1 1 1 0 0 0 1 0 0 0 0 0 0 
+0 0 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 
+0 0 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 
+0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 
+0 0 0 0 0 1 1 1 1 1 0 0 0 0 0 1 0 0 0 0 0 0 
+0 0 0 1 1 1 1 1 1 1 1 0 0 0 1 1 0 0 0 0 1 0 
+0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 
+0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 
+0 1 1 1 1 1 1 1 1 1 1 1 1 0 0 1 0 0 0 0 0 0 
+0 0 1 1 1 1 1 1 1 1 1 0 0 0 0 0 1 0 0 0 0 0 
+0 0 1 1 1 1 1 1 1 1 1 0 0 0 0 0 1 0 0 0 0 0 
+0 0 0 1 1 1 1 1 1 1 0 0 0 0 1 1 0 0 0 0 0 0 
+0 0 0 0 0 0 0 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 1 0 0 
+0 0 0 1 1 1 1 1 1 1 1 1 0 0 1 1 0 0 0 0 1 0 
+0 1 1 1 1 1 1 1 1 1 1 1 0 0 1 1 0 0 0 0 0 0 
+0 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 
+0 0 1 1 1 1 1 1 1 1 1 0 0 0 0 0 1 0 0 0 0 0 
+0 0 1 1 1 1 1 1 1 1 0 0 0 0 0 0 1 0 0 0 0 0 
+0 0 0 0 0 1 1 1 1 1 0 0 1 1 1 1 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 1 1 1 1 1 0 0 0 0 0 0 0 1 1 0 0 0 
+0 0 0 0 1 1 1 1 1 1 0 0 0 0 0 1 1 1 1 0 0 0 
+0 0 0 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 1 0 0 0 0 0 0 
+0 0 0 0 0 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 1 1 1 0 0 1 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 
+0 0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 

HDVF/examples/HDVF/data/data_cubical/Eight_10.pgm

@@ -0,0 +1,87 @@
+P2
+4 12 7
+255
+0 0 0 0 
+0 0 0 0 
+0 0 0 0 
+0 0 0 0 
+0 0 0 0 
+0 0 0 0 
+0 0 0 0 
+0 0 0 0 
+0 0 0 0 
+0 0 0 0 
+0 0 0 0 
+0 0 0 0 
+0 0 0 0 
+0 0 0 0 
+0 1 1 0 
+0 1 1 0 
+0 0 1 0 
+0 0 0 0 
+0 0 0 0 
+0 1 1 0 
+0 1 1 0 
+0 1 1 0 
+0 0 0 0 
+0 0 0 0 
+0 0 0 0 
+0 1 1 0 
+0 1 1 0 
+0 0 0 0 
+0 1 1 0 
+0 1 1 0 
+0 1 1 0 
+0 1 1 0 
+0 0 0 0 
+0 1 1 0 
+0 1 1 0 
+0 0 0 0 
+0 0 0 0 
+0 1 1 0 
+0 0 0 0 
+0 0 0 0 
+0 0 0 0 
+0 1 1 0 
+0 1 1 0 
+0 0 0 0 
+0 0 0 0 
+0 0 0 0 
+0 1 1 0 
+0 0 0 0 
+0 0 0 0 
+0 1 1 0 
+0 1 1 0 
+0 0 0 0 
+0 1 1 0 
+0 1 1 0 
+0 1 1 0 
+0 1 1 0 
+0 0 0 0 
+0 1 1 0 
+0 1 1 0 
+0 0 0 0 
+0 0 0 0 
+0 0 0 0 
+0 0 1 0 
+0 1 1 0 
+0 0 0 0 
+0 0 0 0 
+0 0 0 0 
+0 0 0 0 
+0 1 1 0 
+0 0 1 0 
+0 0 0 0 
+0 0 0 0 
+0 0 0 0 
+0 0 0 0 
+0 0 0 0 
+0 0 0 0 
+0 0 0 0 
+0 0 0 0 
+0 0 0 0 
+0 0 0 0 
+0 0 0 0 
+0 0 0 0 
+0 0 0 0 
+0 0 0 0 

HDVF/examples/HDVF/data/data_cubical/Example_cub_complexes/cub_complex1.cub

@@ -0,0 +1,11 @@
+2
+3 3
+0 0
+1 0
+2 0
+0 1
+2 1
+0 2
+1 2
+2 2
+

HDVF/examples/HDVF/data/data_cubical/Example_cub_complexes/cub_complex2.cub

@@ -0,0 +1,12 @@
+2
+3 3
+0 0
+1 0
+2 0
+0 1
+1 1
+2 1
+0 2
+1 2
+2 2
+

HDVF/examples/HDVF/data/data_cubical/Example_cub_complexes/cub_complex3.cub

@@ -0,0 +1,16 @@
+2
+5 3
+0 0
+1 0
+2 0
+3 0
+4 0
+0 1
+2 1
+4 1
+0 2
+1 2
+2 2
+3 2
+4 2
+