Opened 9 years ago
Closed 2 years ago
#12449 closed enhancement (worksforme)
Improve the way that sage evaluates symbolic functions on basic types
Reported by: | bober | Owned by: | bober |
---|---|---|---|
Priority: | major | Milestone: | sage-duplicate/invalid/wontfix |
Component: | symbolics | Keywords: | gamma function |
Cc: | Snark, eviatarbach | Merged in: | |
Authors: | Jonathan Bober, Burcin Erocal | Reviewers: | Burcin Erocal, Jonathan Bober |
Report Upstream: | N/A | Work issues: | |
Branch: | Commit: | ||
Dependencies: | #4498, #12507, #9130 | Stopgaps: |
Description (last modified by )
Currently when a symbolic function (e.g. gamma(), exp(), sin()) gets a python float, it passes it off to ginac, and let's ginac evaluate it. This can be very slow, e.g.:
sage: timeit('exp(6.0r)') 625 loops, best of 3: 476 µs per loop
This also leads to possibly undesirable inconsistencies across different platforms: for example, ginac ends up calling the local libc tgammal to evaluate the gamma function, and even when this tgammal is evaluated by the same version of eglibc on two different platforms the answer varies depending on whether 80 bit doubles are available.
To fix this, we make Sage check its types which correspond most closely to python type for the function that should be called, and use the functions for those types if available.
Now we have
sage: timeit('exp(6.0r)') 625 loops, best of 3: 1.02 µs per loop
For the specific case of the gamma function, the is also a change here to make RDF use python's math.gamma(). It is more accurate than gsl. It is perhaps not as accurate as eglibc's gamma(), but it should give reliable results on different platforms.
I think this will fix 3 of the 4 failing numerical accuracy doctests on ARM, though the test for binomial(.5r, 5) had to be weakened slightly. (see http://groups.google.com/group/sage-devel/browse_thread/thread/3c8c61ea113ea60c ) binomial() is currently not computed in a very good way, though, so it is reasonable to weaken this test temporarily. (see http://trac.sagemath.org/sage_trac/ticket/12448 )
Finally, I'll remark that these improvements could probably be much better. The above timings should be compared to
sage: timeit('math.exp(6.0r)') 625 loops, best of 3: 112 ns per loop
Apply:
Attachments (1)
Change History (39)
comment:1 Changed 9 years ago by
- Cc Snark added
- Component changed from basic arithmetic to symbolics
- Description modified (diff)
- Status changed from new to needs_review
- Summary changed from Improve the gamma function to Improve the way that sage evaluates symbolic functions on basic types
comment:2 follow-up: ↓ 3 Changed 9 years ago by
- Reviewers set to Burcin Erocal
- Status changed from needs_review to needs_work
comment:3 in reply to: ↑ 2 ; follow-up: ↓ 4 Changed 9 years ago by
Replying to burcin:
AFAICT (without applying the patch and testing), attachment:12449.patch removes the
try/except
block aroundreturn getattr(args[0], self._name)()
, disabling the call toBuiltinFunction.__call__
later on in the function body.
Sorry. I was looking at the code after applying attachment:trac_1173-move_call.patch:ticket:1173. There is no try/except
block in the original code.
Regarding the changes to sage/rings/real_double.pyx
: IIRC, python's math module just calls the underlying libc functions. On ARM won't this introduce new errors?
comment:4 in reply to: ↑ 3 Changed 9 years ago by
Replying to burcin:
Replying to burcin:
AFAICT (without applying the patch and testing), attachment:12449.patch removes the
try/except
block aroundreturn getattr(args[0], self._name)()
, disabling the call toBuiltinFunction.__call__
later on in the function body.Sorry. I was looking at the code after applying attachment:trac_1173-move_call.patch:ticket:1173. There is no
try/except
block in the original code.Regarding the changes to
sage/rings/real_double.pyx
: IIRC, python's math module just calls the underlying libc functions. On ARM won't this introduce new errors?
For the gamma function, at least, python has its own implementation. Of course, that function itself might call some underlying libc functions, so maybe it won't be entirely consistent, now that I think about it.
In the bit of testing that I did, math.gamma was less accurate than tgammal on my computer, but much more accurate than gsl's gamma function.
comment:5 Changed 9 years ago by
- Description modified (diff)
- Status changed from needs_work to needs_review
It turns out that Python implements the gamma function itself:
http://hg.python.org/cpython/file/58bd6a58365d/Modules/mathmodule.c#l226
I uploaded a new patch that changes py_tgamma()
in sage/symbolic/pynac.pyx
to use Python's gamma implementation.
I give a positive review to attachment:12449.patch. If someone can review my patch we can switch this to positive review.
BTW, this will fix #9162. We should close that as a duplicate.
comment:6 Changed 9 years ago by
I uploaded a new version of my patch since the previous function had a typo that make doctests fail.
comment:7 Changed 9 years ago by
- Reviewers changed from Burcin Erocal to Burcin Erocal, Jonathan Bober
- Status changed from needs_review to positive_review
Well, I suppose that I will just go ahead and give it a positive review, since your changes are simple and orthogonal to mine (and I checked that all tests still pass on my machine). I would like to get some confirmation that this fixes things on ARM, but I suppose that we'll find that out eventually; that's a secondary point as far as this patch goes, though, so it isn't that important.
comment:8 Changed 9 years ago by
- Dependencies set to #4498, #12507, #9130
- Status changed from positive_review to needs_work
comment:9 Changed 9 years ago by
- Description modified (diff)
- Status changed from needs_work to needs_review
I've attached a patch rebased to 5.0-beta5. I had no problem applying the patches, so I just made sure all tests still pass.
Maybe I should just go ahead and put it back to positive review, but I'll give Burcin a day or so to take a look at it.
comment:10 Changed 9 years ago by
- Status changed from needs_review to positive_review
Looks good to me.
comment:11 Changed 9 years ago by
- Status changed from positive_review to needs_work
Rebasing to sage-5.0.beta5 isn't sufficient, since #9130 was only merged in sage-5.0.beta6. So, this patch still needs an extra rebase.
comment:12 Changed 9 years ago by
Oops. I guess that explains why I didn't actually need to do anything to rebase. I'll try to remember to update this later today.
comment:13 Changed 9 years ago by
- Description modified (diff)
comment:14 Changed 9 years ago by
- Status changed from needs_work to positive_review
Ok, now it is rebased to 5.0beta7.
comment:15 Changed 9 years ago by
- Status changed from positive_review to needs_work
There is numerical noise. On hawk (OpenSolaris? 06.2009-32):
sage -t --long -force_lib devel/sage/sage/symbolic/function.pyx Warning: divide by zero encountered in divide ********************************************************************** File "/export/home/buildbot/build/sage/hawk-1/hawk_full/build/sage-5.0.beta8/devel/sage-main/sage/symbolic/function.pyx", line 185: sage: cot(complex(1,2)) Expected: (0.03279775553375259-0.9843292264581908j) Got: (0.03279775553375261-0.9843292264581911j) ********************************************************************** sage -t --long -force_lib devel/sage/sage/functions/hyperbolic.py Warning: divide by zero encountered in divide Warning: divide by zero encountered in divide ********************************************************************** File "/export/home/buildbot/build/sage/hawk-1/hawk_full/build/sage-5.0.beta8/devel/sage-main/sage/functions/hyperbolic.py", line 543: sage: float(arccoth(2)) Expected: 0.5493061443340548 Got: 0.5493061443340549 **********************************************************************
comment:16 Changed 9 years ago by
On bsd (OS X 10.6 64-bit):
sage -t --long -force_lib devel/sage/sage/functions/hyperbolic.py Warning: divide by zero encountered in divide Warning: divide by zero encountered in divide ********************************************************************** File "/Users/buildbot/build/sage/bsd-1/bsd_64_full/build/sage-5.0.beta8/devel/sage-main/sage/functions/hyperbolic.py", line 543: sage: float(arccoth(2)) Expected: 0.5493061443340548 Got: 0.5493061443340549 **********************************************************************
comment:17 follow-up: ↓ 18 Changed 9 years ago by
On iras (ia64):
sage -t --long -force_lib devel/sage/sage/rings/arith.py ********************************************************************** File "/home/buildbot/build/sage/cleo-1/cleo_full/build/sage-5.0.beta8/devel/sage-main/sage/rings/arith.py", line 3066: sage: abs(binomial(0.5r, 5) - 0.02734375) < 1e-17r Expected: True Got: False **********************************************************************
comment:18 in reply to: ↑ 17 Changed 9 years ago by
Replying to jdemeyer:
On iras (ia64):
sage -t --long -force_lib devel/sage/sage/rings/arith.py ********************************************************************** File "/home/buildbot/build/sage/cleo-1/cleo_full/build/sage-5.0.beta8/devel/sage-main/sage/rings/arith.py", line 3066: sage: abs(binomial(0.5r, 5) - 0.02734375) < 1e-17r Expected: True Got: False **********************************************************************
Three remarks:
- that test is also one that is failing on ARM ; and there's a trac ticket about binomial's implementation being pretty wrong (#12448).
- should that test use "#tol abs" instead of doing the computation itself?
- shouldn't it be even a "#tol rel" ?
comment:19 Changed 9 years ago by
On cicero (Linux i386):
sage -t --long -force_lib devel/sage/sage/symbolic/function.pyx Warning: divide by zero encountered in divide ********************************************************************** File "/home/buildbot/build/sage/cicero-1/cicero_full/build/sage-5.0.beta8/devel/sage-main/sage/symbolic/function.pyx", line 185: sage: cot(complex(1,2)) Expected: (0.03279775553375259-0.9843292264581908j) Got: (0.032797755533752596-0.9843292264581911j) **********************************************************************
comment:20 Changed 9 years ago by
- Description modified (diff)
- Status changed from needs_work to needs_review
comment:21 Changed 9 years ago by
Ok, here is a version that should quiet those numerical errors. I'm leaving this as need review for the moment instead of putting it back to positive review in case someone wants to revisit the issue of whether or not we might want to make sure that we don't return different results on different machines.
comment:22 Changed 9 years ago by
I think it's better to use RR
(= RealField(53)
) or CC
(=ComplexField(53)
) instead of RDF
or CDF
. That should give more consistent results.
comment:23 follow-up: ↓ 24 Changed 9 years ago by
Missing "-":
sage: abs(float(arccoth(2)) - 0.5493061443340548) < 1e16
comment:24 in reply to: ↑ 23 Changed 9 years ago by
Replying to dsm:
Missing "-":
sage: abs(float(arccoth(2)) - 0.5493061443340548) < 1e16
Oops. That is fixed, and the 16 is changed to a 15, since 16 is actually too large, I think. I also took the opportunity to write a better commit message.
Of course, I'm still not sure if this is the right way to do things.
Changed 9 years ago by
comment:25 follow-up: ↓ 26 Changed 9 years ago by
- Status changed from needs_review to positive_review
Ok, after a brief discussion at SD 40.5 (after Jeroen went to sleep and wasn't there to complain) I'm restoring this to positive review. A few people seemed to agree that if a function is called on a machine floating point type then it is reasonable for it to give machine dependent answers.
I hope that all the tests actually pass now. Also, I've twice updated the patch to remove trailing whitespace, so I hope that I actually got it all this time.
(Of course, Jeroen will still have an opportunity to object...)
comment:26 in reply to: ↑ 25 Changed 8 years ago by
- Status changed from positive_review to needs_work
Replying to bober:
A few people seemed to agree that if a function is called on a machine floating point type then it is reasonable for it to give machine dependent answers.
Obviously the following comment doesn't apply then:
# It is more accurate than gsl's gamma function and # should provide consistant results across platforms.
Also: you should add tests for the infinity stuff that you added.
comment:27 Changed 8 years ago by
* ping *
The issues I mentioned in the last comment are fairly trivial.
comment:28 Changed 7 years ago by
- Milestone changed from sage-5.11 to sage-5.12
comment:29 Changed 7 years ago by
- Cc eviatarbach added
comment:30 Changed 7 years ago by
- Milestone changed from sage-6.1 to sage-6.2
comment:31 Changed 7 years ago by
- Milestone changed from sage-6.2 to sage-6.3
comment:32 Changed 6 years ago by
- Milestone changed from sage-6.3 to sage-6.4
comment:33 Changed 6 years ago by
The patch actually contains code that switches FP sage_tgamma1
from GSL to python_gamma
which might be a good thing, but requires a separate ticket. The timings are identical, so I suspect that the implementation might be too.
comment:34 follow-up: ↓ 36 Changed 6 years ago by
@Jeroen: Is the patch still relevant after the changes you made to BuiltinFunction.__call__
?
comment:35 Changed 6 years ago by
- Status changed from needs_work to needs_info
comment:36 in reply to: ↑ 34 Changed 6 years ago by
comment:37 Changed 3 years ago by
- Milestone changed from sage-6.4 to sage-duplicate/invalid/wontfix
- Status changed from needs_info to needs_review
The original issue seems resolved:
sage: %timeit('exp(6.0r)') 100000000 loops, best of 3: 7.31 ns per loop sage: %timeit('math.exp(6.0r)') 100000000 loops, best of 3: 7.28 ns per loop
comment:38 Changed 2 years ago by
- Resolution set to worksforme
- Status changed from needs_review to closed
Thanks for the patch. I agree that considering float and complex types separately in the fast evaluation path makes sense for speed. However, this causes inconsistencies with the results you get when a value is substituted into the expression.
We still need to modify the
py_tgamma()
function insage/symbolic/pynac.pyx
to fix this:AFAICT (without applying the patch and testing), attachment:12449.patch removes the
try/except
block aroundreturn getattr(args[0], self._name)()
, disabling the call toBuiltinFunction.__call__
later on in the function body.This would be a good occasion to merge my patch from attachment:trac_1173-move_call.patch:ticket:1173. I'll try to rebase your patch on top of that.