test containment of ideals in class MPolynomialIdeal

[sagetrac-mariah]

There seems to be no way to test containment of ideals in the class MPolynomialIdeal in sage.rings.polynomial.multi_polynomial_ideal. One might expect the comparison operators (e.g. I<J ) to do this, but in fact what they do is to compare the Groebner bases as sequences of polynomials, which is counterintuitive. For example:

sage: R.<x,y> = PolynomialRing(QQ)
sage: I=(x*y)*R; J=(x,y)*R; I<J
False
sage: I=(y+1)*R; J=(x,y)*R; I<J
True

This is implemented in the __cmp__ method, which is not up to the task of doing subset comparison, since __cmp__ is only suitable for total orderings.

To do it right would seem to require implementing Python's "rich comparison" methods, __lt__, __gt__, etc.

For example:

from sage.rings.polynomial.multi_polynomial_ideal import MPolynomialIdeal

def IsSubset(I,J):
  for g in I.gens()
    if not g in J: return False
  return True

def IsSuperset(I,J):
  return IsSubset(J,I)

def IsProperSubset(I,J):
  return I!=J and IsSubset(I,J)

def IsProperSuperset(I,J):
  return I!J and IsSuperset(I,J)

setattr(MPolynomialIdeal,'__le__',IsSubset)
setattr(MPolynomialIdeal,'__lt__',IsProperSubset)
setattr(MPolynomialIdeal,'__ge__',IsSuperset)
setattr(MPolynomialIdeal,'__gt__',IsProperSuperset)

With these we now get the expected behavior:

sage: R.<x,y> = PolynomialRing(QQ)
sage: I=(x*y)*R; J=(x,y)*R; I<J
True
sage: I=(y+1)*R; J=(x,y)*R; I<J
False

The patch supplied gives a solution via Groebner bases, and also fixes #12839.

Apply:

1. attachment: trac_12802_no_whitespace.patch

CC: @malb @nthiery

Component: commutative algebra

Keywords: sd40.5, groebner bases, ideals

Author: John Perry

Reviewer: Andrey Novoseltsev, Simon King

Merged: sage-5.4.rc0

Issue created by migration from https://trac.sagemath.org/ticket/12802

[johnperry-math]

Changed keywords from none to sd40.5, groebner bases, ideals

[johnperry-math]

Author: John Perry

[johnperry-math]

comment:1

Attachment: trac_12802_and_12803.patch.gz

The attached file gives a better check for correctness than even the proposed solution, using Groebner bases rather than generators alone. The examples given by mariah produce a correct result on my machine, though for a doctest I use a slightly different solution which her proposed fix wouldn't detect.

This patch also resolve #12839, as far as I can tell.

[novoselt]

Work Issues: documentation

[novoselt]

comment:2

Functions are missing standard documentation blocks.

[johnperry-math]

comment:3

Replying to @novoselt:

Functions are missing standard documentation blocks.

It's even worse than that. For some reason, I'm now finding regressions that I didn't find earlier. I'm working on an updated patch.

[johnperry-math]

comment:4

Okay, so the problems were real, but were discovered because there's a very weird test of elements of quadratic number fields: basically, if a and b are distinct elements of a quadratic number field, then a>b and b>a are always true:

sage: K.<a> = NumberField(x^2 + 1)
sage: a > a + 1
True

That strikes me as really, really odd behavior.

Anyway, I'm about to upload a new patch that undoes the stupidity I did to sage.rings.ideal.py, adds some documentation to the __cmp__ function in Ideal_generic, then takes care of the things that need taking care of down in sage.rings.polynomial.multi_polynomial_idea.py. I've also added documentation. (Sorry about that -- it didn't occur to me because the stuff I was rewriting wasn't documented, either.) It also removes a boatload of whitespace.

[johnperry-math]

Reviewer: novoselt

[johnperry-math]

comment:6

To help the reviewer(s), here are the non-whitespace changes.

1. In sage/rings/ideal.py, I added a docstring & doctest to Ideal_generic's __cmp__.
2. In sage/rings/polynomial/multipolynomial_ideal.py, I
   a. redefined __cmp__ to rely on __lt__, __gt__, and __eq__;
   b. changed __eq__ to test whether each ideal is a subset of the other;
   c. changed __gt__ to test whether other is a subset of self; and
   d. moved the old code for __eq__ to __lt__, where it now checks whether each polynomial in the Groebner basis of self reduces to 0 modulo the Groebner basis of other. Here, I also had to move a try...except loop around, because of a strange AttributeError related to caching, raised only by the doctests of pbori.pyx.

On my machine all the doctests of sage/rings pass with this patch. I haven't doctested everything, though.

[novoselt]

comment:7

Would be better to have whitespaces separately... Anyway: "Decides whether two ideals are equal. The test is performed on the generators." does not sound right, this method does not decide whether ideals are equal. How about "Compare generators of two ideals." instead?

[novoselt]

Changed reviewer from novoselt to Andrey Novoseltsev

[johnperry-math]

comment:8

Replying to @novoselt:

Would be better to have whitespaces separately...

I can do that, if you think it best. I'd rather not, now that I've removed them, but I can.

Anyway: "Decides whether two ideals are equal. The test is performed on the generators." does not sound right, this method does not decide whether ideals are equal. How about "Compare generators of two ideals." instead?

Done.

[novoselt]

comment:9

Apply trac_12802_and_12839.patch

[novoselt]

comment:10

I don't insist on separating whitespace patch now, it is just inconvenient to dig for changes in 100kb. But please make subsequent changes in a new little patch.

1. __cmp__ does not return the result of self == other contrary to what is written. I also think that the test should call cmp directly, since == does not always lead to cmp.
2. In __lt__ the documentation says "We use ideal membership to test whether the generators of each ideal appears in the other." while only one of the containment directions should be checked. Despite of the description of the algorithm in the documentation, the actual code does something else in the beginning and can return numbers or result of cmp instead of True/False.
3. What will __cmp__ of ideals return if neither of them is contained in another? Do we even need it for ideals if there are lt/gt?
4. Documentation of __eq__ is wrong.

[johnperry-math]

comment:11

Replying to @novoselt:

I don't insist on separating whitespace patch now, it is just inconvenient to dig for changes in 100kb. But please make subsequent changes in a new little patch.

Okay, new patch coming. It'll be a bit.

1. __cmp__ does not return the result of self == other contrary to what is written. I also think that the test should call cmp directly, since == does not always lead to cmp.

Got it.

2. In __lt__ the documentation says "We use ideal membership to test whether the generators of each ideal appears in the other." while only one of the containment directions should be checked. Despite of the description of the algorithm in the documentation, the actual code does something else in the beginning and can return numbers or result of cmp instead of True/False.

This was originally code for __cmp__, and when I first started working on it, I had completely forgotten that __cmp__ doesn't work the same as these other special functions. That was one reason for the regression I encountered earlier.

3. What will __cmp__ of ideals return if neither of them is contained in another? Do we even need it for ideals if there are lt/gt?

When I first worked on this, the __cmp__ of Ideal_generic was trumping everything else. That's why I redefined __cmp__, and it's also why I originally did some of this in Ideal_generic.

4. Documentation of __eq__ is wrong.

It took me a while to find what you meant, but I finally saw it.

[johnperry-math]

comment:12

The last line in this if statement that you pointed out:

        if R is not S: # rings are unique
            if type(R) == type(S) and (R.base_ring() == S.base_ring()) and (R.ngens() == S.ngens()):
                other = other.change_ring(R)
            else:
                return cmp((type(R), R.base_ring(), R.ngens()), (type(S), S.base_ring(), S.ngens()))

strikes me as wrong. If the base ring or the number of generators are different, they can't be the same ideal. I'm changing that to return False.

The type is another story. If someone has defined a ring with singular, and wants to compare to a ring defined with cocoa or macaulay, that ought to be possible. Unfortunately, I can't test it: neither Macaulay nor CoCoA build on my Sage 5.0. Are you able to build either?

[novoselt]

comment:13

I think comparison of rings should be done by converting all to a single system where they will be unique. And in general what's the point in ring uniqueness if they may not be unique? If you are not satisfied with R is S check, use R == S and let the meaning of this be handled by ring code, not the ideal comparison.

Note also that there are many failing doctests, the worst ones are with infinite recursion. If it will be necessary to change output in doctests of toric varieties, please do it on top of #13023 which is already positively reviewed.

[johnperry-math]

comment:51

Okay, a diff performed on the results of different patches showed the problems.

Most tests passed before I interrupted the previous run; I've checked that the previous failures now pass. I'm running tests now on the files that were not checked before.

[johnperry-math]

comment:52

All tests pass on my machine.

[simon-king-jena]

comment:53

In order to tell the patch bot what to do:

Apply trac_12802_no_whitespace.patch

[johnperry-math]

comment:54

I just noticed from the patchbot that I had uploaded an older version of the patch, on account of forgetting to sync. If you're testing this, please download the revision. Sorry about that.

[jdemeyer]

Changed work issues from Patch commit messages to Patch commit message

[jdemeyer]

comment:55

The new patch still needs a good commit message. This patch isn't about "no white space", it's about ideals...

[johnperry-math]

Attachment: trac_12802_no_whitespace.patch.gz

should incorporate other patches, without the removal of whitespace

[johnperry-math]

comment:56

Informative commit message added.

[jdemeyer]

Changed work issues from Patch commit message to none

[nbruin]

comment:58

How much precedent do we have to let < and > mean a containment partial ordering? Doesn't the presence of '<' make that "sort" doesn't throw an error and instead produces arbitrary results?

[johnperry-math]

comment:59

Replying to @nbruin:

How much precedent do we have to let < and > mean a containment partial ordering? Doesn't the presence of '<' make that "sort" doesn't throw an error and instead produces arbitrary results?

**This problem predated the patch! ** As the original report points out, < returned a result before, but the result was wrong. There was already a __cmp__ method, so the choices were either to remove it or make sure that something else trumped it. Hence the introduction of __lt__, etc.

As long as there is a __cmp__, there will be an ordering -- and Python supplies a __cmp__ when none is given.

(Or am I wrong?)

[nbruin]

comment:60

Not all python types admit comparison

sage: complex(1,2)<complex(1,3)
TypeError: no ordering relation is defined for complex numbers

but you're right that by default it seems they do:

sage: class B(object):
....:     pass
....: 
sage: b1=B()
sage: b2=B()
sage: b1 < b2
False
sage: b2 < b1
True

That changes with Python 3, though.

[johnperry-math]

comment:61

Apparently, complex is a built-in class; I can't find the source that makes it avoid comparison. However, it makes sense to raise an error in the case where the ideals are incomparable. I wouldn't raise a TypeError, though, but a ValueError. Does that sound acceptable? If so, I'll change the status to Needs Work (AGAIN) and upload a new patch...

[novoselt]

comment:62

What exactly does it mean "incomparable ideals"? If I < J means "Is I a subideal of J?", it always can be answered: if I leaves in a different ring without coercions - the answer is just "no". So I don't think any errors have to be thrown.

I also think that it is useful to be able to sort any mix of anything, so sort should never though an error, even if the result depends on the position of objects in memory.

[nbruin]

comment:63

Replying to @novoselt:

I also think that it is useful to be able to sort any mix of anything, so sort should never though an error, even if the result depends on the position of objects in memory.

That just leads to silently incorrect code. If you want to sort on id you can do sort(L,key=id)

sage: sorted([-10,-11,123,124,1000],key=id)
[123, -10, -11, 124, 1000]
sage: sorted([-10,-11,123,124,1000],key=id)
[124, 1000, 123, -10, -11]
sage: sorted([-10,-11,123,124,1000],key=id)
[-11, 124, 123, -10, 1000]

I'm curious to see a useful application of this.

This is all a bit off-topic, though. There IS precedent in python to use < for partial ordering.

sage: set([2])<set([1,2,3])
True
sage: set([1,2,3])<set([3])
False
sage: set([1,2,3])<set([4,5])
False
sage: set([4,5])<set([1,2,3])
False

so supporting this for ideals should be fine. I'd throw an error if the ideals are "incompatible" for comparison (according to coercion rules for common covering structure, I'd say), but False would be consistent too (it's a partial ordering, after all).

[jdemeyer]

Merged: sage-5.4.rc0