disallow dangerous coercions to RIF

[mstreng]

Remove the following coercions because they make interval arithmetic much less trustworthy.

• RealField --> RealIntervalField
• float --> ComplexIntervalField
• SR --> ComplexIntervalField

For example, one could easily do the following by accident.

    sage: iv = 1 + 2^(-55) + 0. + RIF(1)
    sage: iv.lower(), iv.upper()
    (2.00000000000000, 2.00000000000000)

See also this topic on sage-devel

Already removed in other tickets:

• ComplexField --> ComplexIntervalField

CC: @videlec @mezzarobba @dkrenn @cheuberg @bhutz @JohnCremona @pjbruin

Component: coercion

Keywords: RIF RR interval coercion

Author: Marc Mezzarobba

Branch/Commit: 3c7644f

Reviewer: Vincent Delecroix

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

[sagetrac-tcoffee]

comment:4

I think this change will cause many more problems than it solves.

In practice: I expect this will break a lot of existing code. I often use the RR --> RIF coercion intentionally, because it makes it more convenient to do exactly what I want. For instance:

• In an interval branch-and-bound code, cell widths in a geometric decomposition are expressed as reciprocal powers of two, which can be efficiently stored and used as EXACT floating-point numbers. Combining these in arithmetic operations with intervals gives the same results as explicitly converting them to RIF, without having to do the explicit conversion.

• Interval arithmetic is a fast way to obtain approximate bounds on floating-point calculations over ranges of parameters, by simply replacing the parameter of interest with an interval. To do this without RR --> RIF, one would have to add explicit conversions of all other parameters appearing in combination with this parameter everywhere it occurs.

In theory: I think the motivation is flawed. The change was proposed to avoid accidental coercion of exact ring values to RIF via potentially inexact RR. But the problem in this case is the first step, not the second, and reflects the need for a direct coercion from the exact ring to RIF. In such a case, the lack of such a coercion is the bug. Whereas each value in RR does have a valid coercion to RIF, and the lack of this coercion would also be a bug.

Realistically, I think anyone who is actually

1. using RIF to do reliable computing and
2. coercing values from other exact rings that don't have direct coercions to RIF and are not expressible exactly in RR,

is going to be sufficiently aware of the potential problem that they will handle/check the conversion explicitly. Better to expect this than to inconvenience many users who use this coercion in far less arcane circumstances.

[mstreng]

comment:5

Replying to @sagetrac-tcoffee:

I think this change will cause many more problems than it solves.

Thanks for the comments. I mentioned them in the sage-devel thread where the original discussion took place. Let's move the discussion there.

[dkrenn]

Branch: u/dkrenn/remove-coerce-RR-RIF

[mezzarobba]

comment:10

Does your branch pass the tests? I tried doing something similar a while ago (see u/mmezzarobba/wip/coerce_RR_RIF), but I had to make additional changes, and, more importantly, I concluded that this should wait for #22029 if we don't want comparisons between floating-point numbers and intervals to give very misleading results.

---

New commits:

b30ebc3

disallow coercions from RR to RIF

[mezzarobba]

Commit: b30ebc3

[dkrenn]

comment:11

Replying to @mezzarobba:

Does your branch pass the tests? I tried doing something similar a while ago (see u/mmezzarobba/wip/coerce_RR_RIF), but I had to make additional changes, and, more importantly, I concluded that this should wait for #22029 if we don't want comparisons between floating-point numbers and intervals to give very misleading results.

I am cleaning up my SageMath-versions and found this partial work. As no branch was attached to this ticket, I've uploaded it. So it will, for sure, need some work.

[jdemeyer]

Replying to @mstreng:

For example, one could easily do the following by accident.

    sage: iv = 1 + 2^(-55) + 0. + RIF(1)
    sage: iv.lower(), iv.upper()
    (2.00000000000000, 2.00000000000000)

What's wrong with this? I honestly don't really see a problem with coercion RR -> RIF.

[videlec]

comment:13

Replying to @jdemeyer:

Replying to @mstreng:

For example, one could easily do the following by accident.

    sage: iv = 1 + 2^(-55) + 0. + RIF(1)
    sage: iv.lower(), iv.upper()
    (2.00000000000000, 2.00000000000000)

What's wrong with this? I honestly don't really see a problem with coercion RR -> RIF.

With RIF you always have guaranteed results, with RR you do not. RIF should be as safe as possible to ensure that you actually have proven results. Just compare

sage: r1 = RIF(1/3) + (2.0).cos()
sage: r2 = RIF(1/3) + RIF(2).cos()
sage: r1.endpoints()
(-0.0828135032138091, -0.0828135032138090)
sage: r2.endpoints()
(-0.0828135032138091, -0.0828135032138089)

For me it is similar to the way coercions are implemented between RealField. It goes from more precision to less precision so that you have no artificial big precision in your result.

In conclusion: if any coercion, it should be the other way around (by taking the center).

[jdemeyer]

comment:14

I think a lot depends on how you model mathematically elements of interval fields. There are two ways to interpret them: either as real numbers with some uncertainty or actually as intervals of real numbers.

In the "interval" model, the coercion RR -> RIF makes sense: you identify a real number with a singleton, which is a special case of an interval. On the other hand, the "choose the middle point map" RIF -> RR is neither natural (the middle point is an arbitrary choice) nor a morphism, which are two requirements for a coercion.

In the "real number with uncertainty" model, you could see RIF -> RR as the forgetful map which keeps the real number but forgets about the uncertainty. In that case, this map is more natural than the opposite.

For me, it is pretty clear that MPFI is really meant to model intervals while arb models real numbers with uncertainty.

[mezzarobba]

comment:15

Replying to @jdemeyer:

For me, it is pretty clear that MPFI is really meant to model intervals

I sort of agree, but still, the key property is that operations on MPFI intervals return over-approximations of the set of possible results (in particular, they round the endpoints of the results in the correct direction), and having a coercion from plain floating-point numbers would effectively break that.

What would make sense to me if you want ”intervals of real numbers” with a coercion from the corresponding real numbers themselves would be a separate parent Intervals(R) explicitly parametrized by the parent R = RR, QQ, SR... where the endpoints live.

[jdemeyer]

comment:16

Replying to @mezzarobba:

What would make sense to me if you want ”intervals of real numbers”

RealIntervalField is exactly that for RR

with a coercion from the corresponding real numbers themselves

So, is this an argument in favor of the coercion RR -> RIF?

[mezzarobba]

comment:17

Replying to @jdemeyer:

Replying to @mezzarobba:

What would make sense to me if you want ”intervals of real numbers”

RealIntervalField is exactly that for RR

No: operations in RIF round the upper bound of their result up and the lower bound down. They don't just defer to operations in RR.

with a coercion from the corresponding real numbers themselves

So, is this an argument in favor of the coercion RR -> RIF?

Certainly not.

[jdemeyer]

comment:18

Replying to @mezzarobba:

Replying to @jdemeyer:

Replying to @mezzarobba:

What would make sense to me if you want ”intervals of real numbers”

RealIntervalField is exactly that for RR

No: operations in RIF round the upper bound of their result up and the lower bound down.

Yes, obviously because that's the right thing to do.

They don't just defer to operations in RR.

Technically not, but that is just an implementation detail. RIF is the set of intervals over RR. Whether it's implemented by operations in RR or by a separate library (MPFI in this case) doesn't really matter.

[dkrenn]

comment:19

Replying to @videlec:

What's wrong with this? I honestly don't really see a problem with coercion RR -> RIF.

With RIF you always have guaranteed results, with RR you do not. RIF should be as safe as possible to ensure that you actually have proven results. [...]
For me it is similar to the way coercions are implemented between RealField. It goes from more precision to less precision so that you have no artificial big precision in your result.

In conclusion: if any coercion, it should be the other way around (by taking the center).

Huge +1. (Having reliable numerical computations is a (the) major application of interval arithmetic.)

[mstreng]

comment:21

Replying to @jdemeyer:

For me, it is pretty clear that MPFI is really meant to model intervals

It is not to me. I found various texts (by MPFI developers) that motivate why MPFI exists and they all start by saying things about "guaranteeing results" and "certified enclosures" for floating point real arithmetic. And if an automatic silent coercion mixes a guaranteed result with a non-guaranteed result, this is bad.

People prove theorems using RealIntervalField as "real numbers that are guaranteed to be in a certain interval", and I would personally trust such theorems much less if they use software that allows automatic accidental coercions from RR to RIF. The reason: a mistake or bug could introduce a (non-guaranteed-correct) element of RR and the coercion will silently add it to your RIF element, which now could have an incorrect error bound.

I did not know that people also use the outward-rounding RIF for "arithmetic with intervals" without caring about guaranteed inclusion. I understand that for them it is convenient that the following works:

sage: I = RIF(pi, 2*pi)
sage: p = RR(pi)
sage: I + p
1.?e1

to get the interval [5, 8]. However, it is hardly any extra work to do the direct thing that the coercion model currently does for you, which is:

sage: I + RIF(p)
1.?e1

And we can also implement additional functions, e.g.

sage: I.translate(p)
sage: I.add(p)

So I am very much in favour of removing automatic coercions between RR and RIF. Perhaps with a DeprecationWarning if possible.

[dkrenn]

comment:22

Replying to @mstreng:

So I am very much in favour of removing automatic coercions between RR and RIF. Perhaps with a DeprecationWarning if possible.

If it is considered a bug, then no deprecation is needed...

[jdemeyer]

comment:24

The problem with this ticket is that various places in Sage use this coercion. In particular QQbar:

**********************************************************************
File "src/sage/rings/number_field/number_field.py", line 1513, in sage.rings.number_field.number_field.NumberField_generic.construction
Failed example:
    K.<a> = NumberField(x^3-5,embedding=0)
Exception raised:
    Traceback (most recent call last):
      File "/usr/local/src/sage-config/local/lib/python2.7/site-packages/sage/doctest/forker.py", line 517, in _run
        self.compile_and_execute(example, compiler, test.globs)
      File "/usr/local/src/sage-config/local/lib/python2.7/site-packages/sage/doctest/forker.py", line 920, in compile_and_execute
        exec(compiled, globs)
      File "<doctest sage.rings.number_field.number_field.NumberField_generic.construction[14]>", line 1, in <module>
        K = NumberField(x**Integer(3)-Integer(5),embedding=Integer(0), names=('a',)); (a,) = K._first_ngens(1)
      File "/usr/local/src/sage-config/local/lib/python2.7/site-packages/sage/rings/number_field/number_field.py", line 529, in NumberField
        return NumberField_version2(polynomial=polynomial, name=name, check=check, embedding=embedding, latex_name=latex_name, assume_disc_small=assume_disc_small, maximize_at_primes=maximize_at_primes, structure=structure)
      File "sage/structure/factory.pyx", line 364, in sage.structure.factory.UniqueFactory.__call__ (build/cythonized/sage/structure/factory.c:1964)
        return self.get_object(version, key, kwds)
      File "sage/structure/factory.pyx", line 403, in sage.structure.factory.UniqueFactory.get_object (build/cythonized/sage/structure/factory.c:2197)
        cache_key = _cache_key(cache_key)
      File "sage/misc/cachefunc.pyx", line 642, in sage.misc.cachefunc.cache_key (build/cythonized/sage/misc/cachefunc.c:3189)
        o = cache_key_unhashable(o)
      File "sage/misc/cachefunc.pyx", line 651, in sage.misc.cachefunc.cache_key_unhashable (build/cythonized/sage/misc/cachefunc.c:3503)
        return tuple(cache_key(item) for item in o)
      File "sage/misc/cachefunc.pyx", line 651, in genexpr (build/cythonized/sage/misc/cachefunc.c:3398)
        return tuple(cache_key(item) for item in o)
      File "sage/misc/cachefunc.pyx", line 642, in sage.misc.cachefunc.cache_key (build/cythonized/sage/misc/cachefunc.c:3189)
        o = cache_key_unhashable(o)
      File "sage/misc/cachefunc.pyx", line 651, in sage.misc.cachefunc.cache_key_unhashable (build/cythonized/sage/misc/cachefunc.c:3503)
        return tuple(cache_key(item) for item in o)
      File "sage/misc/cachefunc.pyx", line 651, in genexpr (build/cythonized/sage/misc/cachefunc.c:3398)
        return tuple(cache_key(item) for item in o)
      File "sage/misc/cachefunc.pyx", line 642, in sage.misc.cachefunc.cache_key (build/cythonized/sage/misc/cachefunc.c:3189)
        o = cache_key_unhashable(o)
      File "sage/misc/cachefunc.pyx", line 653, in sage.misc.cachefunc.cache_key_unhashable (build/cythonized/sage/misc/cachefunc.c:3558)
        k = o._cache_key()
      File "sage/structure/element.pyx", line 1059, in sage.structure.element.Element._cache_key (build/cythonized/sage/structure/element.c:9948)
        return(self.parent(),str(self))
      File "sage/structure/sage_object.pyx", line 237, in sage.structure.sage_object.SageObject.__repr__ (build/cythonized/sage/structure/sage_object.c:2778)
        result = repr_func()
      File "sage/rings/real_lazy.pyx", line 748, in sage.rings.real_lazy.LazyFieldElement._repr_ (build/cythonized/sage/rings/real_lazy.c:9768)
        return str(self.approx())
      File "sage/rings/real_lazy.pyx", line 772, in sage.rings.real_lazy.LazyFieldElement.approx (build/cythonized/sage/rings/real_lazy.c:9863)
        return self.eval(self._parent.interval_field())
      File "sage/rings/real_lazy.pyx", line 1640, in sage.rings.real_lazy.LazyAlgebraic.eval (build/cythonized/sage/rings/real_lazy.c:18796)
        roots = self._poly.roots(ring = AA if isinstance(self._parent, RealLazyField_class) else QQbar)
      File "sage/rings/polynomial/polynomial_element.pyx", line 7540, in sage.rings.polynomial.polynomial_element.Polynomial.roots (build/cythonized/sage/rings/polynomial/polynomial_element.c:68834)
        rts = real_roots(self, retval='algebraic_real')
      File "sage/rings/polynomial/real_roots.pyx", line 4155, in sage.rings.polynomial.real_roots.real_roots (build/cythonized/sage/rings/polynomial/real_roots.c:49563)
        return [(AA.polynomial_root(r[1], r[0]), r[2]) for r in intv_roots]
      File "/usr/local/src/sage-config/local/lib/python2.7/site-packages/sage/rings/qqbar.py", line 997, in polynomial_root
        return AlgebraicReal(ANRoot(poly, interval, multiplicity))
      File "/usr/local/src/sage-config/local/lib/python2.7/site-packages/sage/rings/qqbar.py", line 5830, in __init__
        self._interval = self.refine_interval(interval, 64)
      File "/usr/local/src/sage-config/local/lib/python2.7/site-packages/sage/rings/qqbar.py", line 5992, in refine_interval
        return self._real_refine_interval(interval, prec)
      File "/usr/local/src/sage-config/local/lib/python2.7/site-packages/sage/rings/qqbar.py", line 6096, in _real_refine_interval
        new_range = uinfo['endpoint'] - uinfo['value'] / slope
      File "sage/structure/element.pyx", line 1369, in sage.structure.element.Element.__sub__ (build/cythonized/sage/structure/element.c:11538)
        return coercion_model.bin_op(left, right, sub)
      File "sage/structure/coerce.pyx", line 1133, in sage.structure.coerce.CoercionModel_cache_maps.bin_op (build/cythonized/sage/structure/coerce.c:10655)
        raise bin_op_exception(op, x, y)
    TypeError: unsupported operand parent(s) for -: 'Real Field with 64 bits of precision and rounding RNDU' and 'Real Interval Field with 64 bits of precision'

[jdemeyer]

comment:25

Replying to @dkrenn:

If it is considered a bug, then no deprecation is needed...

It's documented behaviour so it's maybe a misfeature but not a bug.

Anyway, it seems that most people think that the coercion RR -> RIF is a misfeature, so I'm willing to remove it in #24371. However, I don't want to deal with the fall-out so if somebody wants to fix [comment:24], please go ahead. It's hard to tell if that's the only issue because that one issue in QQbar causes so many doctest failures.

[mezzarobba]

comment:35

I could not test everything locally due to version issues with pari and giac, so the patchbot may still uncover problems, but this is starting to take shape, and comments are already welcome.

@bhutz: You may want to check the commit touching binary_form_reduce. I don't think I broke anything, but the way intervals are used in the parts I had to modify looks a bit strange to me.

@JohnCremona, @pjbruin: Same remark regarding the changes to elliptic_curves.height.

[JohnCremona]

comment:36

I see no problems with the minor code changes to elliptic curve heights, though I don't think that any of the affected functions was written by me.

[pjbruin]

comment:37

I don't think I contributed much to this precise code either, but the changes certainly look fine to me.

[mezzarobba]

comment:38

Thank you both for your comments. John: indeed, git blame did not tell the whole story. If I understand correctly, you committed the code but it was originally written by Robert Bradshaw.

[JohnCremona]

comment:39

Replying to @mezzarobba:

Thank you both for your comments. John: indeed, git blame did not tell the whole story. If I understand correctly, you committed the code but it was originally written by Robert Bradshaw.

That agrees with my memory. heegner.py was from his thesis, I think, and he rewrote the unpleasant complex interval stuff from a script I had in (I think) gp.

[sagetrac-git]

Changed commit from 88e0764 to 5dcf52c

[sagetrac-git]

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

5c463c8	#15114 remove coercions from RR, RDF, float to RIF
ccf7e2e	#15114 remove coercions from SR, float to CIF
c307736	#15114 force choice of coercion maps from UCF to CC, CIF, ...
97a48a2	#15114 binary_form_reduce: avoid mixing intervals and FP numbers
02308be	#15114 fix doctest tolerance parsing after removal of RR → RIF coercion
7e0d6d9	#15114 minor fixes to elliptic_curves code and tests mixing RR with RIF
5dcf52c	#15114 misc small fixes after removal of dangerous coercions

[bhutz]

comment:41

The changes in binary_form_reduce.py look fine to me. The use of .lower() and .upper() I expect should have been there in the first place and we've just never come across an example where being that careful with the errors mattered.

[mezzarobba]

comment:42

Replying to @bhutz:

The changes in binary_form_reduce.py look fine to me. The use of .lower() and .upper() I expect should have been there in the first place and we've just never come across an example where being that careful with the errors mattered.

Thank you.

All tests pass now, reviewers welcome!

[videlec]

comment:43

This is great. I am positively surprised for seeing only few changes for making this happen!

The only trouble I see is that it might break user code without deprecation notice.

[mezzarobba]

comment:44

Replying to @videlec:

This is great. I am positively surprised for seeing only few changes for making this happen!

The switch to Python 3 (and, I think, also the change to make i an element of ℚ[i] that you recently reviewed) went a long way toward solving the worst issues I found the first time I tried. But it is good news in any case!

The only trouble I see is that it might break user code without deprecation notice.

I agree. But I don't know what to to to ease the transition. Do you see a way to display a deprecation warning when a coercion map is used while keeping the corresponding conversion map working?

That being said, between the issues it indirectly causes and the inconsistency with other interval fields, I would almost call the existence of this coercion a bug.

[videlec]

comment:45

Replying to @mezzarobba:

Replying to @videlec:

This is great. I am positively surprised for seeing only few changes for making this happen!

The switch to Python 3 (and, I think, also the change to make i an element of ℚ[i] that you recently reviewed) went a long way toward solving the worst issues I found the first time I tried. But it is good news in any case!

The only trouble I see is that it might break user code without deprecation notice.

I agree. But I don't know what to to to ease the transition. Do you see a way to display a deprecation warning when a coercion map is used while keeping the corresponding conversion map working?

One possibility would be to support the operations (with a warning) without having the coercions. One cumbersome way to do this is to implement custom __add__, __mul__, etc

class Interval:
    def __add__(self, other):
        if self is an interval and other a floating point:
            warn("don't do that")
            return old_behavior(self, other)
        else:
            return Element.__add__(self, other)

But I don't see a simple way to do it in order to also handle functorial constructions such as vectors or polynomials. Though, given the changes you had to do in the doctests, I think it makes sense to only consider scalars.

That being said, between the issues it indirectly causes and the inconsistency with other interval fields, I would almost call the existence of this coercion a bug.

Agreed, but https://xkcd.com/1172/

[mezzarobba]

comment:46

Replying to @videlec:

One possibility would be to support the operations (with a warning) without having the coercions. One cumbersome way to do this is to implement custom __add__

Unless we do it for floating-point types too (which would be a bit invasive), it will work for interval + float but float + interval will still fail.

[sagetrac-git]

Changed commit from 5dcf52c to ddab847

[sagetrac-git]

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

315ddc2	#15114 remove coercions from RR, RDF, float to RIF
3cc998a	#15114 remove coercions from SR, float to CIF
a51cfe0	#15114 force choice of coercion maps from UCF to CC, CIF, ...
ccd98a7	#15114 binary_form_reduce: avoid mixing intervals and FP numbers
ca84412	#15114 fix doctest tolerance parsing after removal of RR → RIF coercion
584d6e4	#15114 minor fixes to elliptic_curves code and tests mixing RR with RIF
e1aabd6	#15114 misc small fixes after removal of dangerous coercions
ddab847	#15114 restore basic functionality, with deprecation warnings

[mezzarobba]

comment:48

Replying to @mezzarobba:

Replying to @videlec:

One possibility would be to support the operations (with a warning) without having the coercions. One cumbersome way to do this is to implement custom __add__

Unless we do it for floating-point types too (which would be a bit invasive), it will work for interval + float but float + interval will still fail.

Ok, for lack of a better option, I have done that—for RealNumbers only. It was... painful.

[sagetrac-git]

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

674b4f7	#15114 remove coercions from RR, RDF, float to RIF
9abccbb	#15114 remove coercions from SR, float to CIF
3be1c0f	#15114 force choice of coercion maps from UCF to CC, CIF, ...
5fce3b6	#15114 binary_form_reduce: avoid mixing intervals and FP numbers
ee5fd27	#15114 fix doctest tolerance parsing after removal of RR → RIF coercion
ba66ce3	#15114 minor fixes to elliptic_curves code and tests mixing RR with RIF
9f47f4c	#15114 misc small fixes after removal of dangerous coercions
3c7644f	#15114 restore basic functionality, with deprecation warnings

[sagetrac-git]

Changed commit from ddab847 to 3c7644f

[mezzarobba]

comment:50

Rebased to fix a merge conflict due to whitespace changes...

[mezzarobba]

comment:51

Vincent, do you think you will finish the review, or should I look for someone willing to take over?

[videlec]

Reviewer: Vincent Delecroix

[mezzarobba]

comment:53

Great, thank you!

[vbraun]

Changed branch from u/mmezzarobba/15114_RR_to_RIF to 3c7644f