Periodic Alpha Shape 3 is not deterministic

[VincentRouvreau]

Issue Details

When parsing a periodic alpha shape 3, CGAL::dispatch_output does not give a deterministic list of objects The bug is in file https://github.com/CGAL/cgal/blob/master/Periodic_3_triangulation_3/include/CGAL/Periodic_3_Delaunay_triangulation_3.h on the line: CGAL::cpp98::random_shuffle (points.begin(), points.end()); Moreover, a spatial_sort is performed just after. If I remove the random_shuffle, then the object parsing is deterministic.

Source Code

https://gudhi.inria.fr/doc/latest/_alpha_complex_2periodic_alpha_complex_3d_persistence_8cpp-example.html

Environment

• Operating system (Windows/Mac/Linux, 32/64 bits): Ubuntu 18.04 x64
• Compiler: gcc 7.3.0
• Release or debug mode: Release
• Specific flags used (if any):
• CGAL version: 4.13.beta1
• Boost version: 1.65.1
• Other libraries versions if used (Eigen, TBB, etc.): Eigen 3.1.0

[sloriot]

I guess we could pass the random number generator of CGAL using CGAL::get_default_random() to random_shuffle(). That way you would simply have to seed it to get the deterministic behavior. Would it solve your problem?

[VincentRouvreau]

If you have a patch, I am ok to try. As I was explaining, removing the random_shuffle is also a solution, but there may be a good reason why it is here.

[mglisse]

Inside spatial_sort, I see:

        boost::rand48 random;
        boost::random_number_generator<boost::rand48, Diff_t> rng(random);
        CGAL::cpp98::random_shuffle(begin,end,rng);

It seems to me that the code in periodic should be consistent with this one (which always uses the same seed). Actually, looking at the implementation in <CGAL/algorithm.h>, it does seem to use the same technique. So I am confused why it would do something different the second time...

Also, as Vincent mentioned, when is_large_point_set, this random_shuffle is redundant since it is immediately followed by spatial_sort. So it should probably move to the else branch.

[MaelRL]

I guess the random_shuffle is there in the hope that the iterative insertion of points from the randomized point set will make it so we reach 1-cover status faster than inserting points from the non-randomized point set.

Maybe @MoniqueTeillaud knows better!

[MaelRL]

I don't see why the call to CGAL::random_shuffle would create any non-determinism: as @mglisse pointed out, it's initialized with a fixed seed, boost::rand48.

I don't have Gudhi, could you please extract some pure CGAL code that reproduces the issue?

[VincentRouvreau]

Please find my code under the following gist https://gist.github.com/VincentRouvreau/54c5e0c2d1f9a6c94b0f15db37319775 I want to maintain a map from Alpha shape vertex handle to my own vertex handles (here represented as int incremented when a new one is found). I do it 5 times with the same point cloud and it creates 5 output files (0.txt, 1.txt, ...) on each execution. You can see the non deterministic output with diff 0.txt 1.txt, diff 1.txt 2.txt

In file https://github.com/CGAL/cgal/blob/master/Periodic_3_triangulation_3/include/CGAL/Periodic_3_Delaunay_triangulation_3.h if you remove the line: CGAL::cpp98::random_shuffle (points.begin(), points.end()); you will see the result will be deterministic.

Hope it will help.

[lrineau]

@VincentRouvreau commented on Oct 5, 2018, 3:46 PM GMT+2:

In file https://github.com/CGAL/cgal/blob/master/Periodic_3_triangulation_3/include/CGAL/Periodic_3_Delaunay_triangulation_3.h if you remove the line: CGAL::cpp98::random_shuffle (points.begin(), points.end()); you will see the result will be deterministic.

That call should have a third argument, being the seeded random generator.

[lrineau]

@lrineau commented on Oct 5, 2018, 3:54 PM GMT+2:

@VincentRouvreau commented on Oct 5, 2018, 3:46 PM GMT+2:

In file https://github.com/CGAL/cgal/blob/master/Periodic_3_triangulation_3/include/CGAL/Periodic_3_Delaunay_triangulation_3.h if you remove the line: CGAL::cpp98::random_shuffle (points.begin(), points.end()); you will see the result will be deterministic.

That call should have a third argument, being the seeded random generator.

I am wrong. We already said that the default random generator was boost::rand48, and its default seed is always the same.

[MaelRL]

I get the same result (non determinism) whether the random_shuffle is there or not.

The problem is that when you compute alpha shapes, you loop over edges or facets. Those iterators are based on cell addresses in the TDS container and are not deterministic. Thus, you might collect a facet from one cell in one run, and in the next run collect it from the opposite cell. This is what you see in the output of the alpha shapes, e.g. in the diff extract below:

[...]
Cell:
26,28,40,12,
Facet:
-9,31,41,
+31,41,9,
Cell
14,23,42,12,
[...]

The only change you get is sometimes a different orientation for edges/facets, though.

It should be possible to make it deterministic by using time stamps to make the compact container deterministic like it is done in Mesh_3, e.g. Mesh_cell_base_3. More cleanly, there could be some kind of vertex/cell_base in T3 with time stamps (probably extracting the code from Mesh_cell_base_3 and friends towards T3).

[MaelRL]

Putting aside a minor issue, there is one problem preventing the use of timestamps, which is the duplication of vertices when switching to 27 sheets (cf https://github.com/CGAL/cgal/blob/master/Periodic_3_triangulation_3/include/CGAL/Periodic_3_triangulation_3.h#L3613). That's a call to tds.create_vertex(canonical_vh), which is a deep copy of the vertex and thus the timestamp is identical for all the periodic copies. I could do

vh = tds.create_vertex();
vh->set_point(canonical_vertex->point())

but then the two vertices have different (potential) other attributes, which is not very periodic.

If you can ensure that your periodic triangulation will always stay in 1 cover (the troublesome part being the removal of the dummy points) and use time stamps, then periodic alpha shapes will be deterministic.

[MoniqueTeillaud]

@VincentRouvreau May I ask you what are the problems that the non-determinism of iterators is creating to you?

[VincentRouvreau]

@MoniqueTeillaud : In the project GUDHI, for the moment, we have only utilities to compute alpha complex 3d persistence (cf. https://gudhi.inria.fr/doc/2.3.0/examples.html - with specific version, like exact, periodic and/or weighted) based on CGAL Alpha shapes 3. What I am working on at the moment is a class that will encapsulate an alpha complex 3d construction (can be fast/exact and/or weighted/not and/or periodic/not). In order to test my code, the only test I found is to build 2 alpha complexes 3d (one is using Alpha Shapes with Exact Comparisons of Alpha Values, and the other with the non exact comparison). Then I can compare my 2 alpha complexes. I am not really happy with these tests as it is more testing CGAL Alpha shapes 3 than my code, but I did not find another reasonable solution. This works fine for every alpha complexes, except for the periodic versions (weighted/not) because it is non-deterministic. After having a quick look to the CGAL code with @mglisse, it seems this is due to a redundant random_shuffle that could be fixed.

[MoniqueTeillaud]

Thank you for the explanations. Your experiments don't seem to agree with Mael's comment above: "I get the same result (non determinism) whether the random_shuffle is there or not." I am surprised that you don't encounter the same problems with non-periodic alpha-shapes. They are based on the same iterators, which are not more deterministic in the standard Triangulation_3, since, as Mael says above, they are based on comparison of addresses, which can change at each run.

[lrineau]

Successfully tested in CGAL-4.14-Ic-27, and then merged.

[lrineau]

Fixed by #3389.