Opened 3 years ago

Closed 3 years ago

#22552 closed defect (fixed)

2 bugs creating a simple 2-point Polytope

Reported by: was Owned by:
Priority: minor Milestone: sage-8.1
Component: geometry Keywords: base ring, polyhedron, days88, IMA coding sprints
Cc: jipilab, moritz, mkoeppe, vdelecroix Merged in:
Authors: Jean-Philippe Labbé Reviewers: Vincent Delecroix
Report Upstream: N/A Work issues:
Branch: 150c7df (Commits) Commit: 150c7df7962e6ddb06583511e907c1cf6606b863
Dependencies: Stopgaps:

Description

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

Change History (36)

comment:1 Changed 3 years ago by novoselt

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

comment:2 Changed 3 years ago by mkoeppe

  • Cc jipilab moritz mkoeppe added

comment:3 Changed 3 years ago by jipilab

  • Branch set to u/jipilab/22552

comment:4 Changed 3 years ago by jipilab

  • Authors set to Jean-Philippe Labbé
  • Cc vdelecroix added
  • Commit set to a098bc11f9066736c002819b869fd28c1aea310c
  • Keywords base ring polyhedron added
  • Status changed from new to needs_review

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


New commits:

1873e32Added explanation on the base ring
a098bc1Added tests from the ticket

comment:5 Changed 3 years ago by git

  • Commit changed from a098bc11f9066736c002819b869fd28c1aea310c to 0ff050f294e1c9112fe5ff686df9f82cf267a61a

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

0ff050fpep8 cleaning

comment:6 follow-up: Changed 3 years ago by 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`::

comment:7 Changed 3 years ago by vdelecroix

Computations over RDF and RR are different!

comment:8 Changed 3 years ago by vdelecroix

An example

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

comment:9 in reply to: ↑ 6 ; follow-up: Changed 3 years ago by 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.

comment:10 Changed 3 years ago by 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?

comment:11 Changed 3 years ago by 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?

comment:12 in reply to: ↑ 9 Changed 3 years ago by 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.

comment:13 Changed 3 years ago by 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

comment:14 Changed 3 years ago by jipilab

Okay, yes, that sounds good!

comment:15 Changed 3 years ago by jipilab

Should I keep the 3 examples showing the difference also?

comment:16 Changed 3 years ago by 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

comment:17 Changed 3 years ago by 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.

Last edited 3 years ago by jipilab (previous) (diff)

comment:18 Changed 3 years ago by jipilab

  • Status changed from needs_review to needs_work

comment:19 Changed 3 years ago by jipilab

  • Dependencies set to #22605

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

comment:20 Changed 3 years ago by vdelecroix

  • Dependencies #22605 deleted
  • Milestone changed from sage-7.6 to sage-8.1

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

comment:21 Changed 3 years ago by git

  • Commit changed from 0ff050f294e1c9112fe5ff686df9f82cf267a61a to 92b9dd0682ca20c611ab659d0bae0b6fe21f54f1

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

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

comment:22 Changed 3 years ago by git

  • Commit changed from 92b9dd0682ca20c611ab659d0bae0b6fe21f54f1 to c5ff31ec57b9afcc5a799c7ea9a1f3c0a0f0a85c

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

c5ff31eFixed failing tests

comment:23 Changed 3 years ago by jipilab

  • Status changed from needs_work to needs_review

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

comment:24 Changed 3 years ago by 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)).

comment:25 Changed 3 years ago by vdelecroix

comment:26 Changed 3 years ago by vdelecroix

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

comment:27 Changed 3 years ago by vdelecroix

  • Keywords days88 IMA coding sprints added
  • Reviewers set to Vincent Delecroix
  • Status changed from needs_review to needs_work

comment:28 Changed 3 years ago by git

  • Commit changed from c5ff31ec57b9afcc5a799c7ea9a1f3c0a0f0a85c to 50cb0c75e937237d4889272dc2da84bce2108b26

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

50cb0c7Clarified precision and moved tests

comment:29 Changed 3 years ago by jipilab

Dear Vincent,

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

comment:30 Changed 3 years ago by jipilab

  • Status changed from needs_work to needs_review

comment:31 follow-up: Changed 3 years ago by 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?

comment:32 in reply to: ↑ 31 Changed 3 years ago by 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.

comment:33 Changed 3 years ago by git

  • Commit changed from 50cb0c75e937237d4889272dc2da84bce2108b26 to ee54e3008b6fa52797892527ed0f0775ee707e9e

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

ee54e30Rephrased and made it into a warning

comment:34 Changed 3 years ago by vdelecroix

  • Branch changed from u/jipilab/22552 to u/vdelecroix/22552
  • Commit changed from ee54e3008b6fa52797892527ed0f0775ee707e9e to 150c7df7962e6ddb06583511e907c1cf6606b863

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:

150c7df22552: disallow base_ring=RealField(n)

comment:35 Changed 3 years ago by jipilab

  • Status changed from needs_review to positive_review

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

comment:36 Changed 3 years ago by vbraun

  • Branch changed from u/vdelecroix/22552 to 150c7df7962e6ddb06583511e907c1cf6606b863
  • Resolution set to fixed
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.