Opened 8 years ago

Closed 16 months 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 bober)

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)

trac12449.patch (9.5 KB) - added by bober 8 years ago.

Download all attachments as: .zip

Change History (39)

comment:1 Changed 8 years ago by bober

  • Authors set to bober
  • 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: Changed 8 years ago by burcin

  • Authors changed from bober to Jonathan Bober
  • Reviewers set to Burcin Erocal
  • Status changed from needs_review to needs_work

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 in sage/symbolic/pynac.pyx to fix this:

sage: gamma(x).subs(x=6.r)          
120.0

AFAICT (without applying the patch and testing), attachment:12449.patch removes the try/except block around return getattr(args[0], self._name)(), disabling the call to BuiltinFunction.__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.

comment:3 in reply to: ↑ 2 ; follow-up: Changed 8 years ago by burcin

Replying to burcin:

AFAICT (without applying the patch and testing), attachment:12449.patch removes the try/except block around return getattr(args[0], self._name)(), disabling the call to BuiltinFunction.__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 8 years ago by bober

Replying to burcin:

Replying to burcin:

AFAICT (without applying the patch and testing), attachment:12449.patch removes the try/except block around return getattr(args[0], self._name)(), disabling the call to BuiltinFunction.__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 8 years ago by burcin

  • Authors changed from Jonathan Bober to Jonathan Bober, Burcin Erocal
  • 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 8 years ago by burcin

I uploaded a new version of my patch since the previous function had a typo that make doctests fail.

comment:7 Changed 8 years ago by bober

  • 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 8 years ago by jdemeyer

  • Dependencies set to #4498, #12507, #9130
  • Status changed from positive_review to needs_work

This should be rebased to #4498 + #12507 + #9130.

comment:9 Changed 8 years ago by bober

  • 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 8 years ago by burcin

  • Status changed from needs_review to positive_review

Looks good to me.

comment:11 Changed 8 years ago by jdemeyer

  • 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 8 years ago by bober

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 8 years ago by bober

  • Description modified (diff)

comment:14 Changed 8 years ago by bober

  • Status changed from needs_work to positive_review

Ok, now it is rebased to 5.0beta7.

comment:15 Changed 8 years ago by jdemeyer

  • 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 8 years ago by jdemeyer

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: Changed 8 years ago by 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
**********************************************************************

comment:18 in reply to: ↑ 17 Changed 8 years ago by Snark

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:

  1. that test is also one that is failing on ARM ; and there's a trac ticket about binomial's implementation being pretty wrong (#12448).
  2. should that test use "#tol abs" instead of doing the computation itself?
  3. shouldn't it be even a "#tol rel" ?

comment:19 Changed 8 years ago by jdemeyer

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 8 years ago by bober

  • Description modified (diff)
  • Status changed from needs_work to needs_review

comment:21 Changed 8 years ago by bober

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 8 years ago by jdemeyer

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: Changed 8 years ago by dsm

Missing "-":

            sage: abs(float(arccoth(2)) - 0.5493061443340548) < 1e16 

comment:24 in reply to: ↑ 23 Changed 8 years ago by bober

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 8 years ago by bober

comment:25 follow-up: Changed 8 years ago by bober

  • 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 jdemeyer

  • 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 7 years ago by jdemeyer

* ping *

The issues I mentioned in the last comment are fairly trivial.

comment:28 Changed 6 years ago by jdemeyer

  • Milestone changed from sage-5.11 to sage-5.12

comment:29 Changed 6 years ago by eviatarbach

  • Cc eviatarbach added

comment:30 Changed 6 years ago by vbraun_spam

  • Milestone changed from sage-6.1 to sage-6.2

comment:31 Changed 6 years ago by vbraun_spam

  • Milestone changed from sage-6.2 to sage-6.3

comment:32 Changed 5 years ago by vbraun_spam

  • Milestone changed from sage-6.3 to sage-6.4

comment:33 Changed 5 years ago by rws

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: Changed 5 years ago by rws

@Jeroen: Is the patch still relevant after the changes you made to BuiltinFunction.__call__?

comment:35 Changed 5 years ago by rws

  • Status changed from needs_work to needs_info

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

Replying to rws:

@Jeroen: Is the patch still relevant after the changes you made to BuiltinFunction.__call__?

No, this overlaps a lot with #17130.

comment:37 Changed 21 months ago by rws

  • 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 16 months ago by embray

  • Resolution set to worksforme
  • Status changed from needs_review to closed
Note: See TracTickets for help on using tickets.