Repository: giotto-ai/giotto-tda (language: python)

Issue:
[BUG]  fit_transform() of a EuclideanCechPersistence() throws error with 3D point cloud
**Bug description**

Incompatible constructor arguments error is thrown when using fit_transform() of a EuclideanCechPersistence object on a 3-dimensional point cloud. The error does not appear when a 2-dimensional point cloud is used.


**Code to reproduce**

`from gtda.homology import EuclideanCechPersistence`
`import numpy as np`
`X = np.array([[[0.0, 0,0],`
`               [1, 3,0],`
`               [3, 4,0],`
`               [5, 1,0],`
`               [2, 0,0]]])`
`eCech = EuclideanCechPersistence()`
`dgm = eCech.fit_transform(X)`

**Expected behavior**
No error should be thrown, instead, a persistence diagram should be generated


**Actual behaviour**

>   **Traceback (most recent call last):**
 File "/media/lab/Shared/genPDcech.py", line 61, in <module>
    **dgm = eCech.fit_transform(X)**
  File "/home/lab/anaconda3/envs/582/lib/python3.6/site-packages/gtda/utils/_docs.py", line 106, in fit_transform_wrapper
    return original_fit_transform(*args, **kwargs)
  File "/home/lab/anaconda3/envs/582/lib/python3.6/site-packages/sklearn/base.py", line 690, in fit_transform
    return self.fit(X, **fit_params).transform(X)
  File "/home/lab/anaconda3/envs/582/lib/python3.6/site-packages/gtda/homology/simplicial.py", line 1442, in transform
    for x in X)
  File "/home/lab/anaconda3/envs/582/lib/python3.6/site-packages/joblib/parallel.py", line 1041, in __call__
    if self.dispatch_one_batch(iterator):
  File "/home/lab/anaconda3/envs/582/lib/python3.6/site-packages/joblib/parallel.py", line 859, in dispatch_one_batch
    self._dispatch(tasks)
  File "/home/lab/anaconda3/envs/582/lib/python3.6/site-packages/joblib/parallel.py", line 777, in _dispatch
    job = self._backend.apply_async(batch, callback=cb)
  File "/home/lab/anaconda3/envs/582/lib/python3.6/site-packages/joblib/_parallel_backends.py", line 208, in apply_async
    result = ImmediateResult(func)
  File "/home/lab/anaconda3/envs/582/lib/python3.6/site-packages/joblib/_parallel_backends.py", line 572, in __init__
    self.results = batch()
  File "/home/lab/anaconda3/envs/582/lib/python3.6/site-packages/joblib/parallel.py", line 263, in __call__
    for func, args, kwargs in self.items]
  File "/home/lab/anaconda3/envs/582/lib/python3.6/site-packages/joblib/parallel.py", line 263, in <listcomp>
    for func, args, kwargs in self.items]
  File "/home/lab/anaconda3/envs/582/lib/python3.6/site-packages/gtda/homology/simplicial.py", line 1357, in _gudhi_diagram
    cech_complex = CechComplex(points=X, max_radius=self.max_edge_length)
  File "/home/lab/anaconda3/envs/582/lib/python3.6/site-packages/gtda/externals/python/cech_complex_interface.py", line 18, in __init__
    self.thisref = Cech_complex_interface(points, max_radius)
**TypeError: __init__(): incompatible constructor arguments. The following argument types are supported:**
    1. gtda.externals.modules.gtda_cech_complex.Cech_complex_interface(points: List[List[float[2]]], max_radius: float)
Invoked with: array([[0., 0., 0.],
       [1., 3., 0.],
       [3., 4., 0.],
       [5., 1., 0.],
       [2., 0., 0.]]), inf

**Versions**
Linux-5.4.0-74-generic-x86_64-with-debian-buster-sid
Python 3.6.12 |Anaconda, Inc.| (default, Sep  8 2020, 23:10:56) 
[GCC 7.3.0]
NumPy 1.19.2
SciPy 1.5.2
Joblib 1.0.1
Scikit-learn 0.23.2
Giotto-tda 0.4.0

**Additional context**
When the following 2-dimensional point cloud is used, it runs without any error.
`X = np.array([[[0.0, 0],
                [1, 3],
                [3, 4],
                [5, 1],
                [2, 0]]])`



Discussion / hints:
@ektas0330 thanks a lot for opening the issue!  I can indeed reproduce and have determined the reason to be the following line in our bindings:
https://github.com/giotto-ai/giotto-tda/blob/aaa9f871b2ef7b10c391f4803a64d47b34804313/gtda/externals/bindings/cech_complex_bindings.cpp#L27
This restricts the point clouds we currently take to be 2-dimensional in the case of Cech.  I am not the author of these particular bindings and do not remember why this restriction is currently in place.  Perhaps it can be lifted now?  In any case, the current restriction should have been documented.

@MonkeyBreaker do you remember about this?
You're right @ulupo, I don't remember why it was hard coded this way.
But looking at the constructor of `Cech_complex`, it should support without any issue N dimensional points.

We could replace :

```cpp
using Point_cloud = std::vector<std::array<double, 2>>; 
```

By

```cpp
using Point_cloud = std::vector<std::vector<double>>; 
```

And it should work. I'll try it on my side
@MonkeyBreaker thank you for the quick reply.  I'll add a tiny regression test once you are done with this.

@ektas0330 we will be releasing `0.5.0` this week or early next week, so hopefully the bug fix will be included in that version.

Fix patch:
diff --git a/azure-pipelines.yml b/azure-pipelines.yml
index 374c32257..db00842bd 100644
--- a/azure-pipelines.yml
+++ b/azure-pipelines.yml
@@ -46,7 +46,7 @@ jobs:
 
   - task: Cache@2
     inputs:
-      key: '"ccache-wheels-v2021.06.24" | $(Agent.OS) | "$(python.version)"'
+      key: '"ccache-wheels-v2021.06.30" | $(Agent.OS) | "$(python.version)"'
       path: $(CCACHE_DIR)
     displayName: ccache
 
@@ -146,7 +146,7 @@ jobs:
 
   - task: Cache@2
     inputs:
-      key: '"ccache-v2021.06.24" | $(Agent.OS) | "$(python.version)"'
+      key: '"ccache-v2021.06.30" | $(Agent.OS) | "$(python.version)"'
       path: $(CCACHE_DIR)
     displayName: ccache
 
diff --git a/gtda/externals/bindings/cech_complex_bindings.cpp b/gtda/externals/bindings/cech_complex_bindings.cpp
index 89a7e7657..fbf5b9019 100644
--- a/gtda/externals/bindings/cech_complex_bindings.cpp
+++ b/gtda/externals/bindings/cech_complex_bindings.cpp
@@ -21,12 +21,12 @@ namespace Gudhi {
 namespace cech_complex {
 
 class Cech_complex_interface {
+ public:
   using Simplex_tree =
       Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_fast_persistence>;
   using Filtration_value = Simplex_tree::Filtration_value;
-  using Point_cloud = std::vector<std::array<double, 2>>;
+  using Point_cloud = std::vector<std::vector<double>>;
 
- public:
   Cech_complex_interface(const Point_cloud& points,
                          Filtration_value max_radius) {
     cech_complex_ =
@@ -57,7 +57,7 @@ PYBIND11_MODULE(gtda_cech_complex, m) {
       Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_fast_persistence>;
   py::class_<Gudhi::cech_complex::Cech_complex_interface>(
       m, "Cech_complex_interface")
-      .def(py::init<std::vector<std::array<double, 2>>,
+      .def(py::init<Gudhi::cech_complex::Cech_complex_interface::Point_cloud,
                     Simplex_tree::Filtration_value>(),
            "points"_a, "max_radius"_a)
       .def("create_simplex_tree",


Test patch:
diff --git a/gtda/homology/tests/test_simplicial.py b/gtda/homology/tests/test_simplicial.py
index a2c7d9017..fbc2baa8c 100644
--- a/gtda/homology/tests/test_simplicial.py
+++ b/gtda/homology/tests/test_simplicial.py
@@ -41,6 +41,13 @@
                       [0.76291108, -0.68774373],
                       [-1.01932365, -0.05793851]]])
 
+# 3d point cloud example -- same as 2d example but z-coord is 0
+X_pc_3d = np.array([[[2., 2.47942554, 0.],
+                    [2.47942554, 2.84147098, 0.],
+                    [2.98935825, 2.79848711, 0.],
+                    [2.79848711, 2.41211849, 0.],
+                    [2.41211849, 1.92484888, 0.]]])
+
 
 def test_vrp_params():
     metric = 'not_defined'
@@ -368,7 +375,7 @@ def test_cp_not_fitted():
                       [0.45954496, 0.45954497, 1.]]])
 
 
-@pytest.mark.parametrize('X', [X_pc, X_pc_list])
+@pytest.mark.parametrize('X', [X_pc, X_pc_list, X_pc_3d])
 def test_cp_transform(X):
     cp = EuclideanCechPersistence()