Opened 5 years ago

Closed 5 years ago

Last modified 5 years ago

#17130 closed defect (fixed)

Fix coercion bugs in symbolic functions

Reported by: jdemeyer Owned by:
Priority: major Milestone: sage-6.4
Component: symbolics Keywords:
Cc: burcin, rws Merged in:
Authors: Jeroen Demeyer Reviewers: Ralf Stephan
Report Upstream: N/A Work issues:
Branch: abab222 (Commits) Commit:
Dependencies: #17131, #17133 Stopgaps:

Description (last modified by jdemeyer)

This uses coercion correctly:

sage: bessel_Y._eval_(RealField(300)(1), 1.0)
-0.781212821300289

However, it seems that __call__() coerces this result back to the first parent, giving false precision:

sage: bessel_Y(RealField(300)(1), 1.0)
-0.781212821300288684511770043172873556613922119140625000000000000000000000000000000000000000

Same issue with functions which are evaluated using Maxima, which does not support arbitrary precision:

sage: R=RealField(300); elliptic_eu(R(1/2), R(1/8))
0.495073732023201484864216581627260893583297729492187500000000000000000000000000000000000000

The gamma_inc() function also mishandles parents:

sage: gamma_inc(float(0), float(1))
AttributeError: type object 'float' has no attribute 'precision'

Apart from this, this branch also removes lots of boilerplate from _eval_ like

if not isinstance(x, Expression) and not isinstance(y, Expression) and \
        (is_inexact(x) or is_inexact(y)):
    x, y = coercion_model.canonical_coercion(x, y)
    return self._evalf_(x, y, s_parent(x))

by wrapping _eval_ inside the new method _evalf_or_eval_ which automatically does this boilerplate.


Possible follow-ups: #10133, #14766, #16587, #17122, #15200

Change History (58)

comment:1 follow-ups: Changed 5 years ago by kcrisman

Just to clarify, you mean that the precision of the less precise input (1.0) is used, but then for some reason this is sent back to the precision of the first one?

My guess is that this line is to blame.

        # we want to convert the result to the original parent if the input
        # is not exact, so we store the parent here
        org_parent = parent_c(args[0])

I don't have time right now to check this out, but I would not be at all surprised. Or it's the super call just above that, I don't myself 100% understand how the inheritance works, and I keep forgetting when eval versus call are called in the functions. You could try other functions with multiple arguments and see if the same problem happens, or try to reverse the precisions and see what happens to test this hypothesis.

comment:2 in reply to: ↑ 1 Changed 5 years ago by jdemeyer

Replying to kcrisman:

Just to clarify, you mean that the precision of the less precise input (1.0) is used, but then for some reason this is sent back to the precision of the first one?

Yes, the less-precise result is converted to the more-precise parent.

comment:3 in reply to: ↑ 1 Changed 5 years ago by jdemeyer

Replying to kcrisman:

        # we want to convert the result to the original parent if the input
        # is not exact, so we store the parent here
        org_parent = parent_c(args[0])

Especially the args[0] is very suspicious indeed.

comment:4 Changed 5 years ago by jdemeyer

The following patch fixes the problem

  • src/sage/symbolic/function.pyx

    diff --git a/src/sage/symbolic/function.pyx b/src/sage/symbolic/function.pyx
    index 408b6da..5beb01d 100644
    a b cdef class BuiltinFunction(Function): 
    914914        res = super(BuiltinFunction, self).__call__(
    915915                        *args, coerce=coerce, hold=hold)
    916916
    917         # we want to convert the result to the original parent if the input
    918         # is not exact, so we store the parent here
    919         org_parent = parent_c(args[0])
    920 
    921         # convert the result back to the original parent previously stored
    922         # otherwise we end up with
    923         #     sage: arctan(RR(1))
    924         #     1/4*pi
    925         # which is surprising, to say the least...
    926         if org_parent is not SR and \
    927                 (org_parent is float or org_parent is complex or \
    928                 (PY_TYPE_CHECK(org_parent, Parent) and \
    929                     not org_parent.is_exact())):
    930             try:
    931                 return org_parent(res)
    932             except (TypeError, ValueError):
    933                 pass
    934 
    935             # conversion to the original parent failed
    936             # we try if it works with the corresponding complex domain
    937             if org_parent is float or org_parent is complex:
    938                 try:
    939                     from sage.rings.complex_double import CDF
    940                     return complex(CDF(res))
    941                 except (TypeError, ValueError):
    942                     pass
    943             elif hasattr(org_parent, 'complex_field'):
    944                 try:
    945                     return org_parent.complex_field()(res)
    946                 except (TypeError, ValueError):
    947                     pass
    948 
    949917        return res
    950918
    951919    cdef _is_registered(self):

but with quite a bit of doctest failures.

comment:5 Changed 5 years ago by jdemeyer

  • Description modified (diff)

comment:6 Changed 5 years ago by jdemeyer

I think a correct result with the wrong parent is better than a wrong result with the correct parent, so I'm inclined to really remove that block of code and fix individual functions instead.

comment:7 follow-up: Changed 5 years ago by kcrisman

I'm not surprised - the problem is that this kind of code was added at various times to handle certain cases where the output was not the same type as the input (e.g., complex output for real/float input). For instance, does the example mentioned in the code

arctan(RR(1)) 

work properly with this patch? Again, I'm sorry I don't have time to do more than read code for five minutes at this point :(

Edit: I see your most recent comment - yes, that would certainly be helpful, though I have a feeling a LOT of functions would have to be fixed... because some probably implicitly rely on this block but aren't thoroughly tested for unusual input/output.

Edit by jdemeyer: Sorry, edit instead of reply

Last edited 5 years ago by jdemeyer (previous) (diff)

comment:8 follow-up: Changed 5 years ago by kcrisman

Perhaps one could take the coercion mutual parents of all args instead for a minimal change? (Would that work?) I think in principle it's better to handle this all at once rather than having things in each individual function - if possible, of course.

comment:9 in reply to: ↑ 7 Changed 5 years ago by jdemeyer

Replying to kcrisman:

arctan(RR(1)) 

work properly with this patch?

This works as it should, many doctest failures involve Python types:

sage: arctan(float(1))
1/4*pi

(I guess this goes through Pynac)

comment:10 in reply to: ↑ 8 Changed 5 years ago by jdemeyer

Replying to kcrisman:

Perhaps one could take the coercion mutual parents of all args instead for a minimal change?

That still wouldn't fix the case where the function is computed to less precision than the inputs.

I think the following would fix most issues:

  • before calling the function, convert all Python types to the corresponding Sage types (float -> RDF and so on)
  • after calling the function, convert back if needed

comment:11 Changed 5 years ago by vbraun

Not all arguments need to be floating point types, there could be integer or string parameters. So you can't indiscriminately coerce.

IMHO we should just add a clause that if output is float and of lower precision then use that lower precision. In an ideal world we would be able to evaluate everything to arbitrary precision, of course.

comment:12 Changed 5 years ago by jdemeyer

  • Description modified (diff)

comment:13 Changed 5 years ago by jdemeyer

  • Authors set to Jeroen Demeyer

comment:14 Changed 5 years ago by jdemeyer

  • Branch set to u/jdemeyer/ticket/17130
  • Created changed from 10/10/14 09:22:33 to 10/10/14 09:22:33
  • Modified changed from 10/10/14 21:16:07 to 10/10/14 21:16:07

comment:15 Changed 5 years ago by jdemeyer

  • Commit set to 6d107297a5d8867f1a1e9996678ffa441adc8a2d
  • Dependencies set to #17131

New commits:

6d10729Simplify numerical evaluation of BuiltinFunctions

comment:16 Changed 5 years ago by jdemeyer

  • Dependencies changed from #17131 to #17131, #17133

comment:17 Changed 5 years ago by git

  • Commit changed from 6d107297a5d8867f1a1e9996678ffa441adc8a2d to b6e1ed44a663f7410fddb2e3e4c134aa3a0ce8cf

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

0001941Improve accuracy of polytopes.regular_polygon()
24020fePartition().to_exp() should return Sage Integers
b6e1ed4Merge remote-tracking branches 'origin/u/jdemeyer/ticket/17131' and 'origin/u/jdemeyer/ticket/17133' into ticket/17130

comment:18 Changed 5 years ago by jdemeyer

  • Status changed from new to needs_review

comment:19 Changed 5 years ago by jdemeyer

  • Summary changed from Coercion after _eval_() in symbolic functions to Fix coercion bugs in symbolic functions

comment:20 Changed 5 years ago by jdemeyer

  • Description modified (diff)

comment:21 follow-up: Changed 5 years ago by kcrisman

  • Cc burcin added

Wow, that first commit is a nice patch bomb. I hesitate to make so much change to how symbolic functions are evaluated without taking a pretty close look, my apologies for not trying to do that immediately. Also, I think there was an actual reason for calling it inexact and not numerical, though I don't remember offhand.

comment:22 in reply to: ↑ 21 Changed 5 years ago by jdemeyer

Replying to kcrisman:

Wow, that first commit is a nice patch bomb.

Would you be more willing to review the patch if it would be split up? I could try to do that but only if I have confirmation that it will increase the chances of this ticket being reviewed.

comment:23 follow-up: Changed 5 years ago by vbraun

I don't think it is overly long, nor is there a good way of splitting it up. Patch looks good to me.

Its somewhat confusing that is_numerical is completely different from GiNaC::numeric. How about is_approximate, is_limited_precision, or is_floating_point (I'd be happy to call padics "floating point" in this case).

Is it possible to replace is_inexact with is_numerical everywhere? The difference is

sage: is_inexact(pi)
True
sage: sin._is_numerical(pi)
False

I'm not super happy with pi being "inexact", though perhaps there is a technical reason for why we need it.

comment:24 follow-up: Changed 5 years ago by kcrisman

  • Deprecation, folks.
  • No, it isn't long, but I at least need to think about the mechanism by which things are evaluated and want to make sure we don't miss any odd cases. _eval_ has been fairly standard for quite some time.
  • I agree that the naming is problematic no matter what name you choose. I think the rings folks had some good reasons for their choices:
        cpdef bint is_exact(self) except -2:
            """
            Test whether the ring is exact.
    
            .. NOTE::
    
                This defaults to true, so even if it does return ``True``
                you have no guarantee (unless the ring has properly
                overloaded this).
    
            OUTPUT:
    
            Return True if elements of this ring are represented exactly, i.e.,
            there is no precision loss when doing arithmetic.
    
            EXAMPLES::
    
                sage: QQ.is_exact()
                True
                sage: ZZ.is_exact()
                True
                sage: Qp(7).is_exact()
                False
                sage: Zp(7, type='capped-abs').is_exact()
                False
            """
            return True
    
    but the symbolic ring might have inexact stuff in it, so
    Definition:     SR.is_exact(self)
    Source:
        cpdef bint is_exact(self) except -2:
            """
            Return False, because there are approximate elements in the
            symbolic ring.
    
            EXAMPLES::
    
                sage: SR.is_exact()
                False
    
            Here is an inexact element.
    
            ::
    
                sage: SR(1.9393)
                1.93930000000000
            """
            return False
    
    does that make sense?

comment:25 in reply to: ↑ 24 ; follow-up: Changed 5 years ago by jdemeyer

Replying to kcrisman:

  • Deprecation, folks.

What should be deprecated? Is there something which used to work and no longer works now and should be deprecated?

comment:26 in reply to: ↑ 25 ; follow-up: Changed 5 years ago by kcrisman

  • Deprecation, folks.

What should be deprecated? Is there something which used to work and no longer works now and should be deprecated?

I just meant if one changed the name completely, that's all.

comment:27 in reply to: ↑ 23 Changed 5 years ago by jdemeyer

Replying to vbraun:

Its somewhat confusing that is_numerical is completely different from GiNaC::numeric. How about is_approximate, is_limited_precision, or is_floating_point (I'd be happy to call padics "floating point" in this case).

I like none of those 3 alternatives, as there are valid cases for exact integers being "numerical": this isn't yet implemented on this ticket, but for number theoretical functions (binomial, factorial) it absolutely makes sense to consider integers "numerical". So the name should not obviously exclude integers.

In fact, this was one of the reasons to make _is_numerical a method of the function and not a global function (possible support of p-adics is indeed another reason).

Last edited 5 years ago by jdemeyer (previous) (diff)

comment:28 in reply to: ↑ 26 ; follow-up: Changed 5 years ago by jdemeyer

Replying to kcrisman:

I just meant if one changed the name completely, that's all.

If you show me a concrete instance where a deprecation is needed, I'll happily add it.

comment:29 in reply to: ↑ 28 Changed 5 years ago by kcrisman

I just meant if one changed the name completely, that's all.

If you show me a concrete instance where a deprecation is needed, I'll happily add it.

So far it isn't! I was being a little preemptive in terms of the is_* guys. No worries.

comment:30 Changed 5 years ago by jdemeyer

  • Description modified (diff)

comment:31 Changed 5 years ago by rws

  • Cc rws added

comment:32 Changed 5 years ago by jdemeyer

  • Description modified (diff)

comment:33 Changed 5 years ago by jdemeyer

  • Description modified (diff)

comment:34 follow-up: Changed 5 years ago by rws

I have implemented #17151 on top of this patch. Could you please comment on if the code in Func_laguerre._evalf_() needed to determine real_parent can be reduced further, maybe by adding it to this patch, too?

comment:35 in reply to: ↑ 34 Changed 5 years ago by jdemeyer

Replying to rws:

I have implemented #17151 on top of this patch. Could you please comment on if the code in Func_laguerre._evalf_() needed to determine real_parent can be reduced further

Most likely it can be reduced further (why not simply take the parent of x?). I also don't like the name real_parent since it can be complex too.

maybe by adding it to this patch, too?

No. People are already complaining that this patch is too big. But a follow-up ticket certainly makes sense.

Last edited 5 years ago by jdemeyer (previous) (diff)

comment:36 follow-up: Changed 5 years ago by kcrisman

Something else to note (very minor) is that if the top of http://trac.sagemath.org/wiki/symbolics/functions needed some slight updating in order to still be useful, it would be good to do that at the time this was resolved.

comment:37 in reply to: ↑ 36 Changed 5 years ago by jdemeyer

Replying to kcrisman:

Something else to note (very minor) is that if the top of http://trac.sagemath.org/wiki/symbolics/functions needed some slight updating in order to still be useful, it would be good to do that at the time this was resolved.

I don't think anything needs updating. But it would be good though to think more generally about what the API should be for symbolic functions. I consider this ticket as a first step in the right direction.

comment:38 Changed 5 years ago by rws

  • Branch changed from u/jdemeyer/ticket/17130 to public/17130

comment:39 Changed 5 years ago by rws

  • Commit changed from b6e1ed44a663f7410fddb2e3e4c134aa3a0ce8cf to 72659895c652ce20ed1373f95bfea378d9cbe546
  • Reviewers set to Ralf Stephan
  • Status changed from needs_review to needs_work

This is very useful for upcoming code of symbolic functions and contains fixes to functions too. I think I can follow what you do in symbolic. So apart from the doctest fails it's positive from me.

make ptestlong:

sage -t --long src/sage/modular/modform_hecketriangle/abstract_space.py  # 12 doctests failed
sage -t --long src/sage/modular/modform_hecketriangle/readme.py  # 7 doctests failed
sage -t --long src/sage/functions/other.py  # 3 doctests failed

New commits:

382f97aMerge branch 'u/jdemeyer/ticket/17130' of trac.sagemath.org:sage into 6.5beta1
726598917130: reviewer's patch: fix typo

comment:40 Changed 5 years ago by jdemeyer

  • Branch changed from public/17130 to u/jdemeyer/ticket/17130
  • Modified changed from 11/25/14 17:02:40 to 11/25/14 17:02:40

comment:41 Changed 5 years ago by jdemeyer

  • Commit changed from 72659895c652ce20ed1373f95bfea378d9cbe546 to a486db2b801f528a1a00b0dadb4c4c522f21654a
  • Status changed from needs_work to needs_review

New commits:

c47dbd4Fix Trac #17328 again in a better way
a486db2Call the factorial() method of an Integer

comment:42 Changed 5 years ago by rws

  • Status changed from needs_review to positive_review

Passes previously failed tests.

comment:43 Changed 5 years ago by jdemeyer

  • Description modified (diff)

comment:44 Changed 5 years ago by rws

Patchbot failed because of [tutorial ] IOError: [Errno 28] No space left on device

Last edited 5 years ago by rws (previous) (diff)

comment:45 Changed 5 years ago by jdemeyer

  • Description modified (diff)

comment:46 follow-up: Changed 5 years ago by vbraun

I get more numerical noise

sage -t --long src/sage/functions/hyperbolic.py
**********************************************************************
File "src/sage/functions/hyperbolic.py", line 546, in sage.functions.hyperbolic.Function_arccoth.__init__
Failed example:
    float(arccoth(2))
Expected:
    0.5493061443340548
Got:
    0.5493061443340549
**********************************************************************
1 item had failures:
   1 of   9 in sage.functions.hyperbolic.Function_arccoth.__init__
    [156 tests, 1 failure, 1.22 s]
sage -t --long src/sage/functions/hypergeometric.py
    [143 tests, 5.23 s]
sage -t --long src/sage/functions/jacobi.py
    [188 tests, 4.02 s]
sage -t --long src/sage/functions/log.py
**********************************************************************
File "src/sage/functions/log.py", line 661, in sage.functions.log.Function_lambert_w._evalf_
Failed example:
    lambert_w(RDF(-1))
Expected:
    -0.3181315052047642 + 1.3372357014306895*I
Got:
    -0.31813150520476413 + 1.3372357014306895*I
**********************************************************************
File "src/sage/functions/log.py", line 663, in sage.functions.log.Function_lambert_w._evalf_
Failed example:
    lambert_w(float(-1))
Expected:
    (-0.3181315052047642+1.3372357014306895j)
Got:
    (-0.31813150520476413+1.3372357014306895j)
**********************************************************************
1 item had failures:
   2 of   9 in sage.functions.log.Function_lambert_w._evalf_
sage -t --long src/sage/functions/special.py
**********************************************************************
File "src/sage/functions/special.py", line 355, in sage.functions.special.MaximaFunction._evalf_
Failed example:
    f._evalf_(1, I, parent=CDF)
Expected:
    0.8483795707591759 - 0.0742924734216079*I
Got:
    0.8483795707591759 - 0.07429247342160791*I
**********************************************************************
1 item had failures:
   1 of  13 in sage.functions.special.MaximaFunction._evalf_
    [116 tests, 1 failure, 2.35 s]

comment:47 Changed 5 years ago by vbraun

  • Status changed from positive_review to needs_work

comment:48 in reply to: ↑ 46 Changed 5 years ago by rws

Replying to vbraun:

I get more numerical noise

Appears to be system dependent, as all tests in functions/ pass here (OpenSUSE 12.3)

comment:49 Changed 5 years ago by git

  • Commit changed from a486db2b801f528a1a00b0dadb4c4c522f21654a to 9d3cbbdbafb12d1da899abafe0137b4a5639de30

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

9d3cbbdFix numerical noise

comment:50 Changed 5 years ago by jdemeyer

  • Status changed from needs_work to needs_review

comment:51 Changed 5 years ago by rws

  • Status changed from needs_review to positive_review

Still passes tests in functions/ here.

comment:52 Changed 5 years ago by vbraun

  • Status changed from positive_review to needs_work
sage -t --long src/sage/functions/log.py
**********************************************************************
File "src/sage/functions/log.py", line 661, in sage.functions.log.Function_lambert_w._evalf_
Failed example:
    lambert_w(RDF(-1))  # abs tol 1e-16
Expected:
    -0.31813150520476413 + 1.3372357014306895*I
Got:
    -0.31813150520476413 + 1.3372357014306893*I
Tolerance exceeded in 1 of 2:
    + 1.3372357014306895 vs + 1.3372357014306893, tolerance 2e-16 > 1e-16
**********************************************************************
File "src/sage/functions/log.py", line 663, in sage.functions.log.Function_lambert_w._evalf_
Failed example:
    lambert_w(float(-1))  # abs tol 1e-16
Expected:
    (-0.31813150520476413+1.3372357014306895j)
Got:
    (-0.31813150520476413+1.3372357014306893j)
Tolerance exceeded in 1 of 2:
    +1.3372357014306895 vs +1.3372357014306893, tolerance 2e-16 > 1e-16
**********************************************************************
1 item had failures:
   2 of   9 in sage.functions.log.Function_lambert_w._evalf_
    [171 tests, 2 failures, 2.51 s]

comment:53 Changed 5 years ago by git

  • Commit changed from 9d3cbbdbafb12d1da899abafe0137b4a5639de30 to abab222e67040518acee766f747c161a7445bc0a

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

abab222Fix more numerical noise

comment:54 Changed 5 years ago by jdemeyer

  • Status changed from needs_work to positive_review

comment:55 Changed 5 years ago by rws

There may be a followup in a later ticket to deal with an inconsistency, if you agree that the following discrepancy between functions with one and two parameters needs addressing:

sage: i = ComplexField(30).0
sage: gamma(i)
-0.15494983 - 0.49801567*I
sage: gamma(1.0*I)
-0.154949828301811 - 0.498015668118356*I
sage: jacobi_am(2., 2)
0.549820282407325
sage: jacobi_am(2., i)
jacobi_am(2.00000000000000, 1.0000000*I)
sage: jacobi_am(2., 1.0*I)
jacobi_am(2.00000000000000, 1.00000000000000*I)

EDIT: Ah, never mind, jacobi_am takes only reals.

Last edited 5 years ago by rws (previous) (diff)

comment:56 Changed 5 years ago by vbraun

  • Branch changed from u/jdemeyer/ticket/17130 to abab222e67040518acee766f747c161a7445bc0a
  • Resolution set to fixed
  • Status changed from positive_review to closed

comment:57 follow-up: Changed 5 years ago by kcrisman

  • Commit abab222e67040518acee766f747c161a7445bc0a deleted

Question for Jeroen and Ralf - any updating needed in the documentation in sage/symbolic/function_factory.py function function as a result of this ticket? I'd be happy to open a followup if need be.

comment:58 in reply to: ↑ 57 Changed 5 years ago by jdemeyer

Replying to kcrisman:

Question for Jeroen and Ralf - any updating needed in the documentation in sage/symbolic/function_factory.py function function as a result of this ticket?

No, this is only about BuiltinFunctions, there are no functional changes to user-defined functions.

Note: See TracTickets for help on using tickets.