giotto-ai / giotto-tda

A high-performance topological machine learning toolbox in Python
https://giotto-ai.github.io/gtda-docs
Other
845 stars 173 forks source link

Unexpected Cech persistence result returns #669

Open Alessiacosmos opened 1 year ago

Alessiacosmos commented 1 year ago

When using EuclideanCechPersistence to calculate persistence homology, the return result looks wrong.

test code:

import numpy as np
from gtda.homology import EuclideanCechPersistence

exp_data = np.array([[0,0],[2,0],[0.5,1],[2,-1],[1,-1.2]])

# calc PH
homology_dimensions = [0,1] # 0d and 1d PH are calculated.
pers = EuclideanCechPersistence(homology_dimensions=homology_dimensions, n_jobs=6)
diag = pers.fit_transform([exp_data])[0]

The return of diag is:

[[0.         0.78102499 0.        ]
 [0.         0.559017   0.        ]
 [0.         0.50990194 0.        ]
 [0.         0.5        0.        ]
 [0.78102499 1.11803399 1.        ]
 [1.11803401 1.25       1.        ]
 [0.90138781 1.01666667 1.        ]
 [1.         1.00778222 1.        ]
 [1.1280514  1.12805142 1.        ]]

The problems are:

  1. The return result of birth-death seems like radius ($r$) rather than epsilon ($\epsilon=2r$), which doesn't keep the same representation of VR-complex.
  2. When the radius reaches to 1.01666667 or 1.11803399, all 1--simplex(es) in Cech filtration (topological space?) looks have disappeared. So why larger radius are traced? Or do I have any misunderstanding of Cech filtration?

The figure with circle radius=1.01666667 is shown as follows. where the orange edges are edges whose length <=the specific radius 1.01666667 image


I also tried the same point set in GUDHI 3.8(C++ version), and the output is as follows, which looks more reasonable. So I guess if my understanding is right, the backend of CechComplex from GUDHI may need to be updated.

1 0.901388 1.01667
1 1 1.00778
0 0 inf
0 0 0.781025
0 0 0.559017
0 0 0.509902
0 0 0.5

with the core codes:

#include <gudhi/Cech_complex.h>
#include <gudhi/Simplex_tree.h>
#include <gudhi/Persistent_cohomology.h>

#include <CGAL/Epeck_d.h>

#include <iostream>
#include <string>
#include <vector>
#include <algorithm>  // for std::sort

typedef CGAL::Epeck_d< CGAL::Dimension_tag<2> > Kernel;
typedef typename Kernel::Point_d Point;
// simplex tree
typedef Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_fast_persistence>   Simplex_tree;
typedef Simplex_tree::Filtration_value   Filtration_value;

int main(){
    // Type definitions

    // cech complex
    using Point_cloud = std::vector<Point>;
    using Cech_complex = Gudhi::cech_complex::Cech_complex<Kernel, Simplex_tree>;
    // PH
    using Field_Zp = Gudhi::persistent_cohomology::Field_Zp;
    using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology<Simplex_tree, Field_Zp >;

    /**create point cloud**/
    Point_cloud points;
    // [[0,0],[2,0],[0.5,1],[2,-1],[1,-1.2]]
    points.emplace_back(0,0); //0
    points.emplace_back(2,0); //1
    points.emplace_back(0.5,1); //2
    points.emplace_back(2,-1); //3
    points.emplace_back(1,-1.2); //4

    /**Init of a Cech complex from points**/
    Filtration_value max_radius = 3.;
    Cech_complex cech_complex_from_points(points, max_radius);

    /**Simplex tree**/
    Simplex_tree stree;
    cech_complex_from_points.create_complex(stree, 2);

    // ----------------------------------------------------------------------------
    // Compute the persistence diagram of the complex
    // ----------------------------------------------------------------------------
    Persistent_cohomology pcoh(stree);
    // initializes the coefficient field for homology
    pcoh.init_coefficients(2);
    pcoh.compute_persistent_cohomology(0);

    return 0;
}