2 bugs creating a simple 2-point Polytope

[was]

Sara Billey (of Univ of Washington) reported these.

~/Sara$ sage-develop
┌────────────────────────────────────────────────────────────────────┐
│ SageMath version 7.6.beta5, Release Date: 2017-02-26               │
│ Type "notebook()" for the browser-based notebook interface.        │
│ Type "help()" for help.                                            │
└────────────────────────────────────────────────────────────────────┘
┏━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┓
┃ Warning: this is a prerelease version, and it may be unstable.     ┃
┗━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┛
sage: P=Polyhedron(vertices =[(8.3319544851638732, 7.0567045956967727), (6.4876921900819049, 4.8435898415984129)])
sage: P
The empty polyhedron in (Real Field with 57 bits of precision)^2
sage: # WRONG!  It should not be empty.  Indeed, look:
sage: P=Polyhedron(vertices =[(8.3319544851638732, 7.0567045956967727), (6.4876921900819049, 4.84358984159841)])
sage: P
A 1-dimensional polyhedron in RDF^2 defined as the convex hull of 2 vertices
sage: # Also, here's a hub traceback for no reason (as a second but maybe related bug):
sage: P=Polyhedron(vertices =[(8.3319544851638732, 7.0567045956967727), (6.4876921900819049, 4.8435898415984129)], base_ring=RealField(40))
sage: P.plot()
---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
<ipython-input-8-d297b4e23e6b> in <module>()
----> 1 P.plot()
 
/projects/sage/sage-dev/local/lib/python2.7/site-packages/sage/geometry/polyhedron/base.pyc in plot(self, point, line, polygon, wireframe, fill, projection_direction, **kwds)
    694                 return polyhedron.projection()
    695
--> 696         projection = project(self)
    697         try:
    698             plot_method = projection.plot
 
/projects/sage/sage-dev/local/lib/python2.7/site-packages/sage/geometry/polyhedron/base.pyc in project(polyhedron)
    692                 return polyhedron.schlegel_projection()
    693             else:
--> 694                 return polyhedron.projection()
    695
    696         projection = project(self)
 
/projects/sage/sage-dev/local/lib/python2.7/site-packages/sage/geometry/polyhedron/base.pyc in projection(self)
   3759         """
   3760         from .plot import Projection
-> 3761         self.projection = Projection(self)
   3762         return self.projection
   3763
 
/projects/sage/sage-dev/local/lib/python2.7/site-packages/sage/geometry/polyhedron/plot.py in __init__(self, polyhedron, proj)
    492
    493         if polyhedron.ambient_dim() == 2:
--> 494             self._init_from_2d(polyhedron)
    495         elif polyhedron.ambient_dim() == 3:
    496             self._init_from_3d(polyhedron)
 
/projects/sage/sage-dev/local/lib/python2.7/site-packages/sage/geometry/polyhedron/plot.py in _init_from_2d(self, polyhedron)
    748         self.dimension = 2
    749         self._init_points(polyhedron)
--> 750         self._init_lines_arrows(polyhedron)
    751         self._init_area_2d(polyhedron)
    752
 
/projects/sage/sage-dev/local/lib/python2.7/site-packages/sage/geometry/polyhedron/plot.py in _init_lines_arrows(self, polyhedron)
    812             if not obj[i].is_vertex(): continue
    813             for j in range(len(obj)):
--> 814                 if polyhedron.vertex_adjacency_matrix()[i,j] == 0: continue
    815                 if i < j and obj[j].is_vertex():
    816                     l = [obj[i].vector(), obj[j].vector()]
 
/projects/sage/sage-dev/src/sage/misc/cachefunc.pyx in sage.misc.cachefunc.CachedMethodCallerNoArgs.__call__ (/projects/sage/sage-dev/src/build/cythonized/sage/misc/cachefunc.c:13453)()
   2399         if self.cache is None:
   2400             f = self.f
-> 2401             self.cache = f(self._instance)
   2402         return self.cache
   2403
 
/projects/sage/sage-dev/local/lib/python2.7/site-packages/sage/geometry/polyhedron/base.pyc in vertex_adjacency_matrix(self)
   1933             (0, 0, 1, 1, 0) A vertex at (3, 0)
   1934         """
-> 1935         return self._vertex_adjacency_matrix()
   1936
   1937     adjacency_matrix = vertex_adjacency_matrix
 
/projects/sage/sage-dev/local/lib/python2.7/site-packages/sage/geometry/polyhedron/base.pyc in _vertex_adjacency_matrix(self)
    318             M[j, i] = 1
    319
--> 320         face_lattice = self.face_lattice()
    321         for face in face_lattice:
    322             Vrep = face.ambient_Vrepresentation()
 
/projects/sage/sage-dev/src/sage/misc/cachefunc.pyx in sage.misc.cachefunc.CachedMethodCallerNoArgs.__call__ (/projects/sage/sage-dev/src/build/cythonized/sage/misc/cachefunc.c:13453)()
   2399         if self.cache is None:
   2400             f = self.f
-> 2401             self.cache = f(self._instance)
   2402         return self.cache
   2403
 
/projects/sage/sage-dev/local/lib/python2.7/site-packages/sage/geometry/polyhedron/base.pyc in face_lattice(self)
   3408         return Hasse_diagram_from_incidences\
   3409             (atoms_incidences, coatoms_incidences,
-> 3410              face_constructor=face_constructor, required_atoms=atoms_vertices)
   3411
   3412     def faces(self, face_dimension):
 
/projects/sage/sage-dev/local/lib/python2.7/site-packages/sage/geometry/hasse_diagram.pyc in Hasse_diagram_from_incidences(atom_to_coatoms, coatom_to_atoms, face_constructor, required_atoms, key, **kwds)
    180     # Make sure that coatoms are in the end in proper order
    181     tail = [faces[atoms, frozenset([coatom])]
--> 182             for coatom, atoms in enumerate(coatom_to_atoms)]
    183     tail.append(faces[A, frozenset()])
    184     new_order = [n for n in new_order if n not in tail] + tail
 
KeyError: (frozenset([]), frozenset([0]))
sage: 

See public worksheet: ​https://cloud.sagemath.com/projects/53b9d6b6-ce2c-4007-843a-257cc01bf65b/files/Sara/Polygon%20Bug.sagews

[novoselt]

If I understood correctly, #18220 tries to fix this and other problems by disallowing inexact fields completely...

[jipilab]

I added tests and some explanation in the documentation of the polyhedron constructor. This complements #18220 also inside the constructor class.

---

New commits:

​1873e32	Added explanation on the base ring
​a098bc1	Added tests from the ticket

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​0ff050f

pep8 cleaning

[vdelecroix]

Do you really want to keep this behavior?

If the input consists of decimal numbers and the `base_ring` is not specified,
the base ring is set to be the `RealField` with the precision given by the
minimal precision of the input. If it has 53 bits of precision, the constructor
converts automatically the base ring to `RDF`::

[vdelecroix]

Computations over RDF and RR are different!

[vdelecroix]

An example

sage: float("2.0") ** -2048r == 0.0r
True
sage: RR(2.) ** -2048r == RR.zero()
False

[jipilab]

Replying to vdelecroix:

Do you really want to keep this behavior?

If the input consists of decimal numbers and the `base_ring` is not specified,
the base ring is set to be the `RealField` with the precision given by the
minimal precision of the input. If it has 53 bits of precision, the constructor
converts automatically the base ring to `RDF`::

No, not really...

But at least describing the current behavior may be useful to understand how the constructor works. Is that reasonable?

Once this behavior is changed, I'd just delete the paragraph.

[jipilab]

This ticket only looks at the bugs, now with #18220 merged, the behavior is changed and the "bugs" have now an explanation.

What else should be done here?

[jipilab]

Should the "Do what I mean" be changed completely?

Or should the base_ring passed by the constructor to the parent *always* be RDF and not RR?

[vdelecroix]

Replying to jipilab:

Replying to vdelecroix:

Do you really want to keep this behavior?

If the input consists of decimal numbers and the `base_ring` is not specified,
the base ring is set to be the `RealField` with the precision given by the
minimal precision of the input. If it has 53 bits of precision, the constructor
converts automatically the base ring to `RDF`::

No, not really...

But at least describing the current behavior may be useful to understand how the constructor works. Is that reasonable?

When documented, a bug becomes a feature! What about

.. WARNING::

    When used with input in the realf field with 53 bits of mantissa precision
    the input is converted without notice into the machine floating point numbers.
    This behavior might cause some undesirable effects.

[vdelecroix]

An example of undesirable feature

sage: a = RR(2.0)**-2048
sage: b = RR(3.0)**-2048
sage: Polyhedron([[a],[b]])
A 0-dimensional polyhedron in RDF^1 defined as the convex hull of 1 vertex

[jipilab]

Okay, yes, that sounds good!

[jipilab]

Should I keep the 3 examples showing the difference also?

[jipilab]

Further, the warning bugs me a bit...

It does not tell the complete truth about the behaviour: you only need 1 entry to have a 53 bit mentissa to get into floating points not necessary all of them... This is even worse I find! This is somehow the source of the weird behavior described in this ticket. The fact that one value had 53 bit of mentissa made it possible to create the object, and I would say lead to the belief that the polyhedron is well defined.

sage: c = 1.00000000000000
sage: d = 1.000000000000001
sage: Polyhedron([[c],[d]])
A 0-dimensional polyhedron in RDF^1 defined as the convex hull of 1 vertex

[jipilab]

I know not everyone is a fan of issuing a warning, but IMHO, I would find it a useful hint to the user if:

-it gave entries with more than 53 bits of mentissa -it also gave entries with exactly 53 bits of mentissa (not in ZZ or QQ)

To say:

Warning: some entries had more precision than 53 bits and got converted to 53 bits of precision.

This would not be issued if no entries had more precision than 53 bits.

[jipilab]

We should probably add a dependancy to #22605 because the latter touches the constructor and floats...

[vdelecroix]

This is indeed fixed in 8.1.beta2 in which #22605 is merged

sage: P=Polyhedron(vertices =[(8.3319544851638732, 7.0567045956967727), (6.4876921900819049, 4.8435898415984129)])
Traceback (most recent call last):
...
ValueError: no appropriate backend for computations with Real Field with 57 bits of precision
sage: P=Polyhedron(vertices =[(8.3319544851638732, 7.0567045956967727), (6.4876921900819049, 4.84358984159841)])
sage: P
A 1-dimensional polyhedron in RDF^2 defined as the convex hull of 2 vertices

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​92b9dd0

Merge branch 'u/jipilab/22552' of git://trac.sagemath.org/sage into 22552

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​c5ff31e

Fixed failing tests

[jipilab]

Tests seem to pass on 8.1.beta2 but it would be nice to have feedback on the added documentation.

[vdelecroix]

You should be more explicit about what is meant by precision

sage: (1.1).precision()
53
sage: (1.143234134123123).precision()
54
sage: (1.14323413412312123).precision()
60

It is not the number of digits but closer to min(53, bits(x)).

[vdelecroix]

This is not the syntax for test blocks

.. TESTS:

see ​http://doc.sagemath.org/html/en/developer/coding_basics.html#documentation-strings

[vdelecroix]

Since there is four examples in the examples, I don't see why the 'TESTS' section is relevant in your commit.

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​50cb0c7

Clarified precision and moved tests

[jipilab]

Dear Vincent,

I made some corrections addressing your comments. Have a look.

[vdelecroix]

This sentence is definitely too complicated

If the input consists of decimal numbers and the `base_ring`
is not specified, the base ring is set to be the `RealField`
with the precision given by the minimal bit precision of the
input and 53 (the precision of `RDF`). Then, if the obtained
minimum was 53 bits of precision, the constructor converts
automatically the base ring to `RDF`. Otherwise, it returns an error

I would try to be more informative

Be careful when you construct polyhedra with floating point
numbers. The only available backend for such computation is
`cdd` which uses machine floating point numbers which have
have limited precision (24 bits on 32 bits architecture and
53 bits on 64 bits architecture). Sage behavior 
... etc ... 

BTW, what is happening on 32 bits?

[jipilab]

BTW, what is happening on 32 bits?

It seems that binary64 format from

​https://en.wikipedia.org/wiki/IEEE_754

has 53 bits of precision no matter if it was compiled in 32 or 64 bits processors. I'm no expert here but it seems to be safe to say that a Double has 53 bits in 32 and 64 bits.

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​ee54e30

Rephrased and made it into a warning

[vdelecroix]

A (bad) undocumented feature

sage: Polyhedron(vertices =[(8.3319544851638732, 7.0567045956967727), (6.4876921900819049, 4.8435898415984129)], base_ring=RR)
A 1-dimensional polyhedron in RDF^2 defined as the convex hull of 2 vertices

As agreed with JP, I wrote a commit that makes the above an error.

---

New commits:

​150c7df

22552: disallow base_ring=RealField(n)

[jipilab]

The last changes look good to me. I give my green light on the last changes.