-
-
Notifications
You must be signed in to change notification settings - Fork 473
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Fast cardinality method for IntegerVectorsModPermutationGroup #36814
Conversation
sage: I = IntegerVectorsModPermutationGroup(G, max_part=1) | ||
sage: I.cardinality() | ||
11 | ||
""" |
This comment was marked as resolved.
This comment was marked as resolved.
Sorry, something went wrong.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Yes, some fancier groups would be nice in TESTS. For example, it is nice to demonstrate and test that PermutationGroup([(1,2),(3,4)])
(independent swaps, so a group of order 4) is indeed different from PermutationGroup([[(1,2),(3,4)]])
(simultaneous swap, so a group of order 2).
For a "long" test, we could take all subgroups of SymmetricGroup(k)
for k=3 (6 groups) and k=4 (30 groups). Perhaps a few direct products of small symmetric and cyclic groups.
For extensive tests we have to use the sgs
parameter to IntegerVectorsModPermutationGroup
, because otherwise the tests will fail for groups that only differ by the domain, due to #36527 which is still unfixed (this PR does not fix it and it is unclear to me how and when it would be fixed).
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I am +1000 on this contribution. I'll provide few tests for "fancy groups", if you don't mind.
m = self._max_part # Max of one entry, -1 for no limit | ||
if m == -1: | ||
m = d # Any entry cannot exceed total | ||
if k == 0: |
This comment was marked as resolved.
This comment was marked as resolved.
Sorry, something went wrong.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Umm... I'd rather not. We are really talking about integers here, I think it is much more natural to compare to zero than to treat the integer as a pseudo Boolean. Is there an established rule for using not k
pro k == 0
when k
is really and conceptually an integer value?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
It is faster and more common (at least within Sage's library) to do this. Most integers end up not being Python int
's as well, meaning coercion gets involved.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Faster my left foot. The comparison is done once (as in, "1 time") per invocation of cardinality
. The difference is what, 40 nanoseconds? Did you measure?
sage: V=IntegerVectorsModPermutationGroup(SymmetricGroup(5),100)
sage: timeit("V.fast_cardinality()", number=1000, precision=5)
1000 loops, best of 3: 2.4201 ms per loop
sage: timeit("V.faster_cardinality()", number=1000, precision=5)
1000 loops, best of 3: 2.4255 ms per loop
sage: timeit("V.cardinality()", number=1, precision=5)
1 loop, best of 3: 10.837 s per loop
Here cardinality
is the old iterative version. fast_cardinality
is the version with k==0
and faster_cardinality
is exactly same code with not k
. Perhaps the not k
version is faster by a few nanoseconds, but I am not seeing it here. We are replacing an old 10-second algorithm with a new 2.4-millisecond algorithm, and are we now seriously talking about shaving nanoseconds?
Since when did SageMath programmers become interested in shaving nanoseconds anyway? Quoting
https://doc.sagemath.org/html/en/faq/faq-contribute.html
Fast code is valuable, but clean, readable code is important. In the mathematics community, inaccurate results are unacceptable. ... Don Knuth observes that: “We should forget about small efficiencies, say about 97% of the time: premature optimization is the root of all evil.”
Replacing the clean, semantically transparent integer comparison k==0
with a pseudo-Boolean comparison not k
is, in my humble opinion, a prime example of a futile attempt at premature optimization that obscures the meaning and works as an open invitation to bugs. And it is against the FAQ that I quoted above.
I ask again: Is this not k
style for integer comparison endorsed in some official documentation, or is it just a habit that SageMath programmers have grown fond of, for better or worse?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
For the record, yes, I have repeatedly measured. Not every coding convention (in fact, many of them are not) is in some official documentation. Not everything is called "1 time" such as the method being called in a tight loop. This is also not premature optimization but a standard micro-optimization; premature optimization is writing very complicated code for something that is not a bottleneck (as one has not profiled code). For example, another oddity of Python is tuple([i for i in range(5)])
is faster than tuple(i for i in range(5))
, so we often use the former even though the latter is more "clean" (in some sense of the word).
I hope you will continue to want to be an active and collaborative part of the SageMath developers community. Yet, I feel I must also point out our code of conduct and that your tone is close to not adhering to those standards.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Well, if that thing is a "standard micro-optimization" routinely applied in SageMath code, I must say I am sorry to hear that.
With all due respect I disagree that it is a good idea to apply it all over the place, or here. I have expressed the reasons for my stance and I will leave it at that.
I will not be applying that kind of thing here, but if someone else wants to, then that's their right.
This comment was marked as resolved.
This comment was marked as resolved.
Sorry, something went wrong.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
@mantepse Thank you for your considerations. Informed opinions, agreeable or otherwise, are always welcome.
Agreeing is nice, just I'm not sure about how much agreement is achievable. I hope we agree that replacing a default algorithm with a
If the not k
kind of optimization is deemed crucial and mandatory in SageMath, well, I am happy to report that I won't be bothered. I have plenty of other things to do in research. How mandatory it is, well, that is up to the community.
This comment was marked as resolved.
This comment was marked as resolved.
Sorry, something went wrong.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
So, could you simply reshape the pull request otherwise as discussed
Of course, because those other things make sense.
Re: variable names, I think my variable names and in-code comments are already more descriptive and less terse than what I have seen in other SageMath code (again I am very sorry to see that). But I am not one to be against code clarity, as you may have noticed.
m = d # Any entry cannot exceed total | ||
if k == 0: | ||
# Special case: Empty vectors cannot have nonzero sum. | ||
if d == 0 or d is None: |
This comment was marked as resolved.
This comment was marked as resolved.
Sorry, something went wrong.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Exactly, better not.
else: | ||
return Integer(0) | ||
|
||
if d is None: |
This comment was marked as resolved.
This comment was marked as resolved.
Sorry, something went wrong.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
No, you have to multiply the power series. To wit,
sage: G=SymmetricGroup(2)
sage: Z=G.cycle_index()
sage: m=2
sage: V=IntegerVectorsModPermutationGroup(G, max_part=m)
sage: V.cardinality()
6
sage: sum(c*prod((m+1)**cycle_len for cycle_len in cycle_type) for cycle_type, c in Z)
9
sage: list(V)
[[0, 0], [1, 0], [2, 0], [1, 1], [2, 1], [2, 2]]
The correct count is 6 as we can see from the list.
This comment was marked as resolved.
This comment was marked as resolved.
Sorry, something went wrong.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Yes, that is correct. This is because the function-counting series is then just a k*m
-degree polynomial, and the sum of its coefficients is the value of the polynomial at x==1
. (There is probably a more direct combinatorial way to see this also.) This adds a bit of complication/duplication to the code, but I guess it is worth it because it reduces to simple integer arithmetic, avoiding power series computations altogether in the d is None
case.
series_prod((1 - x**((m+1)*cyclen)) / (1 - x**cyclen) | ||
for cyclen in cyctype) | ||
* coeff | ||
for (cyctype, coeff) in Z) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I suggest to do for cycle_type, c in G.cycle_index())
(and do cycle_len
and cycle_type
instead of the other abbreviations.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I do not see much difference in clarity between cyctype
and cycle_type
in this context, and I deliberately chose the shorter one for the benefit of the reader, to better differentiate a local variable from the cycle_type
method of PermutationGroupElement
.
But it may be argued that a well-equipped reader is well aware of what is a local variable, so yes, we can use a longer nonabbreviated name. What I don't understand is why in your other code snippet you are then asking for the opposite: for the descriptive coeff
to be replaced with the terse c
?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I don't mind much. I read cyclen
as a weird spelling of the plural of cycle
(I'm Austrian, i.e. German speaking), and I was confused (even though I know my Pólya-Redfield quite well). Since separating words with underscores is very common, I suggested it. It is really as minor as it gets.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Oh, I did not see that way of reading :-)
I have absolutely no issues against using full words and underscore here, I will change it the way you like. But is it ok to keep coeff
? You do not want it to be tersified into c
?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I don't mind. Use what fits you.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Oh good. I was afraid that the single letter was another house rule directed against code clarity. (You know, the compiler may be 2 nanoseconds faster that way, and you never know when someone is going to call the compiler a billion times in a tight loop.)
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Please try to keep this a happy, fun and productive space. Please resist making jokes about other opinions.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I'm not joking (though I must admit that recently I have sometimes had difficulty keeping a straight face). I have seen so much terse SageMath code with single-letter variable names and nary a comment, that it did cross my mind that it is a house rule.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I believe we can all agree that there are many (subjective) descriptions of what constitutes sufficient code clarity. Because of that, we can work together as a community without judgement to improve the overall utility of the software.
@@ -775,6 +784,79 @@ def __iter__(self): | |||
else: | |||
return SF.elements_of_depth_iterator(self._sum) | |||
|
|||
def cardinality(self): | |||
r"""Return the number of integer vectors in the set. |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This should be
r"""
Return ...
(note the newline) per sage conventions
Few examples counting 0-1 vectors, e.g.
Indeed, PGL(2,q), is 3-transitive:
but not 4-transitive. The stabiliser of a 3-set is |
I separated the For the benefit of those who care about actual speed, this is what the line profiler says: With If the performance of |
Hei @jukkakohonen - I've sent you an invite to a sagemath triage team - in practice this means less red tape here on GitHub in sagemath repos, not really any extra obligations. |
Nice extra examples @dimpase ! Indeed 0-1 vectors are an important use case. (Graphs in general are 0-1 vectors, when viewed that way.) Do you mean these as "examples" or "tests"? |
Thanks, I have been wondering how many times someone else just has to press the "go ahead" button every time I push a new version and they need to be linted.... |
In some special cases there are faster ways of finding the cardinality than full-fledged Pólya counting. I'm wondering how many special cases would be sensible here. Empty domain is already handled specially.
As an example for 1., consider the trivial group over 10 elements, and let the sum argument be 20. The general method takes roughly 600 microseconds, but calculating the cardinality by the binomial coefficient As an example for 2., consider the following setting where only 2 balls are distributed over 10 boxes, with some nontrivial symmetry.
Even here Pólya counting is not much worse than iterating. With the same group and sum=4, cardinality becomes 207 and iterating is already slower than Pólya. So it would really be very small cardinalities that are useful to calculate by iteration. I guess that 1. is a very common use case (if you construct permutation groups e.g. from the automorphisms of some big graphs, it often happens that the group is in fact trivial). It would be easy to recognize and handle it as a special case here, so that the user does not have to. Furthermore one could try to recognize some special cases of groups (like cyclic or symmetric) and invoke special code, but I somehow think it would be an overkill here. If some particular use case involves lots of groups of some specific structure, allowing perhaps faster counting, then probably the user is in a better position to know about the structure and to exploit it if needed. |
# Special case: Empty vectors cannot have nonzero sum. | ||
if d == 0 or d is None: | ||
return Integer(1) | ||
else: |
This comment was marked as resolved.
This comment was marked as resolved.
Sorry, something went wrong.
if self._sum is not None and self._sum > 0: | ||
# No empty vector can have positive sum. | ||
return iter(()) | ||
else: |
This comment was marked as resolved.
This comment was marked as resolved.
Sorry, something went wrong.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Not sure what you aim to achieve by these "else-unnecessary" comments. You probably know I am not going to change the code to something I consider inferior, in terms of style, clarity and reliability.
This comment was marked as resolved.
This comment was marked as resolved.
Sorry, something went wrong.
Thanks, I added the PGL example. |
if d is not None and self.has_trivial_group() and m >= d: | ||
# Simple calculation with stars and bars. | ||
if d is not None and m >= d and all(g.is_one() for g in G.gens()): | ||
# Trivial group (should use G.is_trivial() as soon as |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
You can and should make the pr for ·is_trivial· a dependency, then we don't need to change things again.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Right, but for that I need to find out how to make a dependency between PR's (never did that before). The current code works fine as it is, regardless of the is_trivial
pr, so I'm not convinced it pays off to complicate things by creating dependencies.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
There is a heading "Dependencies" at the very top of this PR, and there you simply write #36873
. Then, locally, you checkout the other PR (eg., gh pr checkout 36873
, switch back to your PR, do git merge perm_gps/is_trivial
, and push the result.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Thanks, I hope I did it right.
@@ -763,18 +788,192 @@ def __iter__(self): | |||
[3, 1, 1, 2] | |||
[3, 0, 2, 2] | |||
[2, 2, 2, 1] | |||
|
|||
Check that :issue:`36681` is fixed:: | |||
sage: I = IntegerVectorsModPermutationGroup( |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Indentation is wrong here: Check ...
needs to go 4 spaces to the left. Also, add a new blank line after ::
Finally, for some editors (eg., mine), the following style is preferrable::
sage: G = PermutationGroup([], domain=[])
sage: I = IntegerVectorsModPermutationGroup(G, sum=0)
(i.e., avoid ....:
)
(same thing below)
Check that :issue:`36681` is fixed:: | ||
sage: I = IntegerVectorsModPermutationGroup( | ||
....: PermutationGroup([], domain=[]), sum=0) | ||
sage: list(iter(I)) # Single empty vector |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Personally, I'd remove the two comments, here and two lines below, because they say exactly what the output is saying anyway.
|
||
EXAMPLES:: | ||
|
||
With a trivial group all vectors are canonical:: |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
As above, the indentation needs to be fixed. Also, in this case, EXAMPLES::
should be EXAMPLES:
(because there will be no indent after). [this is sphinx, not me :-)]
from sage.rings.rational_field import QQ | ||
from sage.rings.integer import Integer | ||
from sage.misc.misc_c import prod | ||
from sage.arith.misc import binomial |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I don't care much, but usually the imports are at the top of the file: https://peps.python.org/pep-0008/#imports
21 | ||
|
||
With two interchangeable elements, the smaller one | ||
ranges from zero to n//2:: |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
n//2
needs double backticks (for the html doc to come out correctly)
Return the number of integer vectors in the set. | ||
|
||
The number is computed using the cycle index theorem, which is | ||
faster than listing the vectors. |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I'd actually remove , which is faster than listing the vectors
, because it carries no useful information, at least not for the user.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I think it can be left in - perhaps moved into ..NOTE
. It does carry useful info (if you only care about the number, don't try to list vectors)
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I don't really mind, but I don't see the point either - it is obvious that you wouldn't implement a formula if enumeration is faster, isn't it?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I don't understand why a user would not benefit from being told what is fast and what is not. Certainly that is more relevant from a user's viewpoint than which particular algorithm is being used. So perhaps we remove the cycle index theorem then.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I disagree, but let's leave it at that, then. Let's not argue about seven words.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I will remove the information as requested.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
For some reason the default rank
and unrank
still talk about speed.
This is the default (brute force) implementation from the category
"EnumeratedSets()" which can be used when the method "__iter__" is
provided. Its complexity is O(r), where r is the rank of "obj".
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Let me just add my personal reason why I find "uses the Redfield–Pólya theorem" valuable information and "faster than plain enumeration" not.
Before that, please note that I don't really mind, because I know the code and the math behind it quite well, so it doesn't affect me in this case, and I see that one can have other opinions. In particular, you can leave it as it is, or modify it, as you wish.
The first gives me some rough information on how the result is obtained. This can be helpful if I have to understand the code (eg., whether I can adapt it to other, similar classes) and also what I can expect from it performance-wise. For example, if the doc string says: "We essentially use linear search to find the element of given rank.", this tells me immediately that, if I know a custom unranking algorithm, I will be able to improve the performance. Also, the code might use some ad-hoc tricks to make it slightly faster, which will almost certainly make it harder to decipher what the underlying algorithm really is.
I don't see any information I could gain from "faster than plain enumeration". If it would be slower than plain enumeration, I wouldn't use a (nontrivial) formula, would I?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I will remove the information as requested.
No, please don't do this. It is absolutely normal, and preferred, for Sage documentation to explain what mathematical tools are behind the implementation.
If a user wants to understand how a piece of code works, it's not a good idea to leave it totally unexplained, just like a blackbox. Sage has users from all walks of mathematical and CS life, I don't know what percentage of them heard about the cycle index theorem, perhaps 5% or less. It's not a topic that comes up in every maths-related undergraduate curriculum.
I'd even like to ask for a reference to the cycle index theorem to be added to the docstring.
It could be just
:wikipedia:`Cycle_index`
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I don't see any information I could gain from "faster than plain enumeration". If it would be slower than plain enumeration, I wouldn't use a (nontrivial) formula, would I?
Just imagine you didn't know what the cycle index theorem is all about. The choice of the algorithm might have been like "we chose an algorithm we knew, and implemented, just because we could". Only after an explanation that it's a huge speedup it'd made sense.
as explained in my comments on the PR
There is a linter failure, other than that looks good to me. |
Thank you for the extra polish! |
Documentation preview for this PR (built with commit 11bbef9; changes) is ready! 🎉 |
Fast cardinality method for IntegerVectorsModPermutationGroup
This patch fixes #36787 by implementing a
cardinality
method forIntegerVectorsModPermutationGroup_with_constraints
. The methodcalculates the cardinality using the the Polya enumeration
theorem (also known as the cycle index theorem).
It is faster than the the default implementation, which iterates through the full set to
find the cardinality.
Incidentally this PR fixes #36681 so that cardinality and iter no longer crash in empty-domain situations.
📝 Checklist
⌛ Dependencies
is_trivial
from that PR