Implementation of Cohen-Macaulay test for simplicial complexes

[Christian Stump]

Simplicial complexes were lacking a method to test Cohen-Macaulayness.

In order to implement this method, some other methods are improved, namely a hash function is implemented, and _enlarge_subcomplex has become faster.

For convinience, I also added a face_iterator.

Remark: the new line

int_facets = set( a.set().intersection(f_set) for a in new_facets ) 

in _enalarge_subcomplex improved speed for computing the homology by 65% (in the example I looked at -- needs to be verified). This method still has the potential to be speeded a lot, and it is responsible for a lot cpu time when computing the homology.

I also added a second version using parallel tests on multiple cpus.

---

Apply only: trac_11523-cohen_macaulay_complex-cs-jhp-review-ts.patch​

[John Palmieri]

I'm sorry, I forgot about this ticket.

Overall, I think this looks pretty good. Some comments and questions:

• I had to rebase this against 4.7.2.alpha2.
• why have both is_cohen_macaulay and is_cohen_macaulay_parallel? When would you ever want to use the first one? With ncpus=1, the parallel one is about as fast as the first, so should we get rid of the first one?
• on the other hand, when I run the parallel version with ncpus=1, I sometimes see an error:

  Exception OSError: (10, 'No child processes') in <generator object __call__ at 0x10d20ee10> ignored
  

  and then the computation proceeds. Should we try to catch this and hide it? Can we do that if it says "ignored"?
• in my testing (timeit('simplicial_complexes.NotIConnectedGraphs(5,2).homology()')), I saw no change in the time required to do homology computations from the old version to the new one, but if you have other ideas for speeding things up, please implement them.
• the method face_iterator has no documentation.

I don't know if this was ready for review, but I'm attaching a new version of the patch which makes some of these changes: it renames is_cohen_macaulay to _is_cohen_macaulay_serial (so it's hidden), renames the parallel version to is_cohen_macaulay. I wonder if we should set ncpus to be 2 by default instead of 1? I'll try testing on a machine with just one cpu to see what effect that has.

[Christian Stump]

Replying to jhpalmieri:

Thanks for looking at it again, I hope we can get it into sage soon, as it was pretty much done. The only reason why I stopped working was that I thought I can still speed things up, but I never found the time to do so. I vote for leaving it as it is for now, and if someone (maybe me) needs a faster version, we get back to it.

Overall, I think this looks pretty good. Some comments and questions:

• why have both is_cohen_macaulay and is_cohen_macaulay_parallel? When would you ever want to use the first one? With ncpus=1, the parallel one is about as fast as the first, so should we get rid of the first one?

that was just a left-over, as I first implemented is_cohen_macaulay and then started to play with the parallel version.

• on the other hand, when I run the parallel version with ncpus=1, I sometimes see an error:

  Exception OSError: (10, 'No child processes') in <generator object __call__ at 0x10d20ee10> ignored
  

  and then the computation proceeds. Should we try to catch this and hide it? Can we do that if it says "ignored"?

If I remember right, this happens not only for 1 cpu, but if a negative answer is returned without waiting the other parallel sessions to be finished. I agree that we should catch that error!

• in my testing (timeit('simplicial_complexes.NotIConnectedGraphs(5,2).homology()')), I saw no change in the time required to do homology computations from the old version to the new one, but if you have other ideas for speeding things up, please implement them.

See my comment at the beginning.

• the method face_iterator has no documentation.

I will add it!

I don't know if this was ready for review, but I'm attaching a new version of the patch which makes some of these changes: it renames is_cohen_macaulay to _is_cohen_macaulay_serial (so it's hidden), renames the parallel version to is_cohen_macaulay. I wonder if we should set ncpus to be 2 by default instead of 1? I'll try testing on a machine with just one cpu to see what effect that has.

Thanks again, I will update your patch as soon as I have done the changes...

Best, Christian

[Christian Stump]

I updated a new version. The only changes are that I import QQ as we need that, and I simplified the code for is_cohen_macaulay slightly. I saw that you added the documentation for face_iterator already, thanks!

[Christian Stump]

Replying to stumpc5:

I updated a new version. The only changes are that I import QQ as we need that, and I simplified the code for is_cohen_macaulay slightly. I saw that you added the documentation for face_iterator already, thanks!

It is rebased against 4.7.2.alpha3.

[John Palmieri]

This is mostly good. The only problems are with the parallel version: if I run with more than one process, I get a message like

sage: X.is_cohen_macaulay(2)
Killing any remaining workers...
False

If I run it with one process, I get

sage: X.is_cohen_macaulay(1)
Killing any remaining workers...
Exception OSError: (10, 'No child processes') in <generator object __call__ at 0x10ca39a00> ignored
False

We should be able to catch the error. Can we also suppress the message "Killing any remaining workers..."?

[Christian Stump]

Replying to jhpalmieri:

This is mostly good. The only problems are with the parallel version: if I run with more than one process, I get a message like

sage: X.is_cohen_macaulay(2)
Killing any remaining workers...
False

If I run it with one process, I get

sage: X.is_cohen_macaulay(1)
Killing any remaining workers...
Exception OSError: (10, 'No child processes') in <generator object __call__ at 0x10ca39a00> ignored
False

We should be able to catch the error. Can we also suppress the message "Killing any remaining workers..."?

I tried to catch it, but it it already caught somewhere inside the decorator, so I don't know how to suppress the output, do you?

Feel free to update the patch if you see a way to do it!

[John Palmieri]

I'm sorry this has languished for so long. Here is a new patch which turns off the message "Killing remaining workers...". I don't know how to get rid of the "No child processes warning", so I mentioned its possible presence in the docstring, so people won't be too concerned if they see it. I also removed the serial version altogether.

[Travis Scrimshaw]

Because of the hashing of simplicial complexes, this is in conflict with #12587. I've started a topic on sage-devel on the mutability of SimplicialComplex.

[Christian Stump]

Replying to tscrim:

Because of the hashing of simplicial complexes, this is in conflict with #12587. I've started a topic on sage-devel on the mutability of SimplicialComplex.

Thanks for bringing this up. I will modify this ticket here as soon as the other ticket is ready.

[Travis Scrimshaw]

Hey Christian,

Just wondering about the status of this patch now that #12587 is (code-wise) finalized.

Thanks,
Travis

[Christian Stump]

Okay, this is now ready for its final review, I guess... Thanks Travis for bringing it up again!

[Travis Scrimshaw]

Hey Christian and John,

This did not apply against #12587, so I did the very minor rebase in my review patch (it contains the entire patch; sorry if I overstepped here). I also made some minor docstring changes: some formatting and moved the note about the exception in a .. NOTE:: block. If you agree with the changes, feel free to set this to positive review.

Thanks for your work on this, I plan on putting this to good use,
Travis

For patchbot:

Apply only: trac_11523-cohen_macaulay_complex-cs-jhp-review-ts.patch

[Christian Stump]

Replying to tscrim:

Thanks for your work on this, I plan on putting this to good use,

great!

The patch looks good to me, so I am marking it positive...

Best, Christian

[Jeroen Demeyer]

Is the work issue still relevant?

[Travis Scrimshaw]

Replying to jhpalmieri:

Thanks for looking at it again, I hope we can get it into sage soon, as it was pretty much done. The only reason why I stopped working was that I thought I can still speed things up, but I never found the time to do so. I vote for leaving it as it is for now, and if someone (maybe me) needs a faster version, we get back to it.

It seems like it's not.