#17130 closed defect (fixed)
Fix coercion bugs in symbolic functions
Reported by:  jdemeyer  Owned by:  

Priority:  major  Milestone:  sage6.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 )
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.
Change History (58)
comment:1 followups: ↓ 2 ↓ 3 Changed 6 years ago by
comment:2 in reply to: ↑ 1 Changed 6 years ago by
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 lessprecise result is converted to the moreprecise parent.
comment:3 in reply to: ↑ 1 Changed 6 years ago by
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 6 years ago by
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): 914 914 res = super(BuiltinFunction, self).__call__( 915 915 *args, coerce=coerce, hold=hold) 916 916 917 # we want to convert the result to the original parent if the input918 # is not exact, so we store the parent here919 org_parent = parent_c(args[0])920 921 # convert the result back to the original parent previously stored922 # otherwise we end up with923 # sage: arctan(RR(1))924 # 1/4*pi925 # 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 pass934 935 # conversion to the original parent failed936 # we try if it works with the corresponding complex domain937 if org_parent is float or org_parent is complex:938 try:939 from sage.rings.complex_double import CDF940 return complex(CDF(res))941 except (TypeError, ValueError):942 pass943 elif hasattr(org_parent, 'complex_field'):944 try:945 return org_parent.complex_field()(res)946 except (TypeError, ValueError):947 pass948 949 917 return res 950 918 951 919 cdef _is_registered(self):
but with quite a bit of doctest failures.
comment:5 Changed 6 years ago by
 Description modified (diff)
comment:6 Changed 6 years ago by
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 followup: ↓ 9 Changed 6 years ago by
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?
That example works, the problem lies with Python types:
sage: arctan(float(1)) 1/4*pi
comment:8 followup: ↓ 10 Changed 6 years ago by
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 6 years ago by
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 6 years ago by
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 6 years ago by
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 6 years ago by
 Description modified (diff)
comment:13 Changed 6 years ago by
comment:14 Changed 6 years ago by
 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 6 years ago by
 Commit set to 6d107297a5d8867f1a1e9996678ffa441adc8a2d
 Dependencies set to #17131
New commits:
6d10729  Simplify numerical evaluation of BuiltinFunctions

comment:16 Changed 6 years ago by
 Dependencies changed from #17131 to #17131, #17133
comment:17 Changed 6 years ago by
 Commit changed from 6d107297a5d8867f1a1e9996678ffa441adc8a2d to b6e1ed44a663f7410fddb2e3e4c134aa3a0ce8cf
Branch pushed to git repo; I updated commit sha1. New commits:
0001941  Improve accuracy of polytopes.regular_polygon()

24020fe  Partition().to_exp() should return Sage Integers

b6e1ed4  Merge remotetracking branches 'origin/u/jdemeyer/ticket/17131' and 'origin/u/jdemeyer/ticket/17133' into ticket/17130

comment:18 Changed 6 years ago by
 Status changed from new to needs_review
comment:19 Changed 6 years ago by
 Summary changed from Coercion after _eval_() in symbolic functions to Fix coercion bugs in symbolic functions
comment:20 Changed 6 years ago by
 Description modified (diff)
comment:21 followup: ↓ 22 Changed 6 years ago by
 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 6 years ago by
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 followup: ↓ 27 Changed 6 years ago by
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 followup: ↓ 25 Changed 6 years ago by
 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='cappedabs').is_exact() False """ return True
but the symbolic ring might have inexact stuff in it, soDefinition: 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 ; followup: ↓ 26 Changed 6 years ago by
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 ; followup: ↓ 28 Changed 6 years ago by
 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 6 years ago by
Replying to vbraun:
Its somewhat confusing that
is_numerical
is completely different fromGiNaC::numeric
. How aboutis_approximate
,is_limited_precision
, oris_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).
comment:28 in reply to: ↑ 26 ; followup: ↓ 29 Changed 6 years ago by
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 6 years ago by
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 6 years ago by
 Description modified (diff)
comment:31 Changed 6 years ago by
 Cc rws added
comment:32 Changed 6 years ago by
 Description modified (diff)
comment:33 Changed 6 years ago by
 Description modified (diff)
comment:34 followup: ↓ 35 Changed 6 years ago by
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 6 years ago by
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 determinereal_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 followup ticket certainly makes sense.
comment:36 followup: ↓ 37 Changed 6 years ago by
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 6 years ago by
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 6 years ago by
 Branch changed from u/jdemeyer/ticket/17130 to public/17130
comment:39 Changed 6 years ago by
 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:
382f97a  Merge branch 'u/jdemeyer/ticket/17130' of trac.sagemath.org:sage into 6.5beta1

7265989  17130: reviewer's patch: fix typo

comment:40 Changed 6 years ago by
 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 6 years ago by
 Commit changed from 72659895c652ce20ed1373f95bfea378d9cbe546 to a486db2b801f528a1a00b0dadb4c4c522f21654a
 Status changed from needs_work to needs_review
comment:42 Changed 6 years ago by
 Status changed from needs_review to positive_review
Passes previously failed tests.
comment:43 Changed 6 years ago by
 Description modified (diff)
comment:44 Changed 6 years ago by
Patchbot failed because of [tutorial ] IOError: [Errno 28] No space left on device
comment:45 Changed 6 years ago by
 Description modified (diff)
comment:46 followup: ↓ 48 Changed 6 years ago by
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 6 years ago by
 Status changed from positive_review to needs_work
comment:48 in reply to: ↑ 46 Changed 6 years ago by
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 6 years ago by
 Commit changed from a486db2b801f528a1a00b0dadb4c4c522f21654a to 9d3cbbdbafb12d1da899abafe0137b4a5639de30
Branch pushed to git repo; I updated commit sha1. New commits:
9d3cbbd  Fix numerical noise

comment:50 Changed 6 years ago by
 Status changed from needs_work to needs_review
comment:51 Changed 6 years ago by
 Status changed from needs_review to positive_review
Still passes tests in functions/
here.
comment:52 Changed 6 years ago by
 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 1e16 Expected: 0.31813150520476413 + 1.3372357014306895*I Got: 0.31813150520476413 + 1.3372357014306893*I Tolerance exceeded in 1 of 2: + 1.3372357014306895 vs + 1.3372357014306893, tolerance 2e16 > 1e16 ********************************************************************** File "src/sage/functions/log.py", line 663, in sage.functions.log.Function_lambert_w._evalf_ Failed example: lambert_w(float(1)) # abs tol 1e16 Expected: (0.31813150520476413+1.3372357014306895j) Got: (0.31813150520476413+1.3372357014306893j) Tolerance exceeded in 1 of 2: +1.3372357014306895 vs +1.3372357014306893, tolerance 2e16 > 1e16 ********************************************************************** 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 6 years ago by
 Commit changed from 9d3cbbdbafb12d1da899abafe0137b4a5639de30 to abab222e67040518acee766f747c161a7445bc0a
Branch pushed to git repo; I updated commit sha1. New commits:
abab222  Fix more numerical noise

comment:54 Changed 6 years ago by
 Status changed from needs_work to positive_review
comment:55 Changed 6 years ago by
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.
comment:56 Changed 6 years ago by
 Branch changed from u/jdemeyer/ticket/17130 to abab222e67040518acee766f747c161a7445bc0a
 Resolution set to fixed
 Status changed from positive_review to closed
comment:57 followup: ↓ 58 Changed 6 years ago by
 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 6 years ago by
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 BuiltinFunction
s, there are no functional changes to userdefined functions.
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.
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 wheneval
versuscall
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.