Opened 10 years ago
Last modified 5 years ago
#9130 closed enhancement
Access to beta function — at Version 24
Reported by: | kcrisman | Owned by: | burcin |
---|---|---|---|
Priority: | major | Milestone: | sage-5.0 |
Component: | symbolics | Keywords: | special function, pynac, sd35.5 Cernay2012 |
Cc: | benjaminfjones | Merged in: | |
Authors: | Karen T. Kohl, Burcin Erocal | Reviewers: | Benjamin Jones, Burcin Erocal |
Report Upstream: | N/A | Work issues: | |
Branch: | Commit: | ||
Dependencies: | Stopgaps: |
Description (last modified by )
Although Maxima has the beta function, Sage doesn't:
sage: a, b, c = var('a b c') sage: assume(a > 0) sage: assume(b > 0) sage: x = var('x') sage: beta_dist = x**(a-1) * (1 - x)**(b-1) sage: c = integral(beta_dist, x, 0, 1) sage: c beta(a, b) sage: c(a=.5,b=.5) beta(0.500000000000000, 0.500000000000000) sage: c(a=.5,b=.5).n() --------------------------------------------------------------------------- TypeError Traceback (most recent call last) /Users/karl-dietercrisman/<ipython console> in <module>() /Users/karl-dietercrisman/Desktop/sage-4.4.2/local/lib/python2.6/site-packages/sage/symbolic/expression.so in sage.symbolic.expression.Expression.n (sage/symbolic/expression.cpp:17042)() TypeError: cannot evaluate symbolic expression numerically
This *is* is Ginac, though, and there is even room for defining it in symbolic/expression.pyx . It probably is also included in some of our other libraries, as a standard special function.
Apply
Change History (33)
comment:1 Changed 10 years ago by
comment:2 Changed 10 years ago by
I wasn't suggesting which of the several packages in Sage should be used for numerical evaluation, though mpmath did come to mind :-)
I don't think that beta(0.5,0.5) would work, given that
sage: gamma(0.5) 1.77245385090552 sage: gamma(1/2) sqrt(pi)
but beta(1/2,1/2) becoming pi should work fine once we have a symbolic wrapper (with or without Ginac):
sage: maxima_console() Maxima 5.21.1 http://maxima.sourceforge.net using Lisp ECL 10.4.1 Distributed under the GNU Public License. See the file COPYING. Dedicated to the memory of William Schelter. The function bug_report() provides bug reporting information. (%i1) beta(1/2,1/2); (%o1) %pi
Please do try to add this! We definitely often get email asking for various special functions both symbolically and numerically. Also, the more examples we have, the easier it is to finish off the rest of them by cut and paste.
comment:3 follow-up: ↓ 4 Changed 10 years ago by
GiNaC has a beta function, so this can probably be solved simply by wrapping that. See #8864 for an example.
Though I don't know why the beta()
function in sage/symbolic/expression.pyx
is commented. Maybe there is something I'm missing.
Fredrik, it would be great if you can do this. I'd be happy to answer questions if anything goes wrong.
comment:4 in reply to: ↑ 3 ; follow-up: ↓ 5 Changed 10 years ago by
Replying to burcin:
GiNaC has a beta function, so this can probably be solved simply by wrapping that. See #8864 for an example.
Though I don't know why the
beta()
function insage/symbolic/expression.pyx
is commented. Maybe there is something I'm missing.
I think the same reason the psi and psi2 ones are commented - when those were implemented, they didn't notice that they had been commented earlier. This was probably pretty early in the conversion, maybe when William was dealing with CLN (whatever that is)?
comment:5 in reply to: ↑ 4 Changed 10 years ago by
Replying to kcrisman:
Replying to burcin:
GiNaC has a beta function, so this can probably be solved simply by wrapping that. See #8864 for an example.
Though I don't know why the
beta()
function insage/symbolic/expression.pyx
is commented. Maybe there is something I'm missing.I think the same reason the psi and psi2 ones are commented - when those were implemented, they didn't notice that they had been commented earlier. This was probably pretty early in the conversion, maybe when William was dealing with CLN (whatever that is)?
psi()
and psi2()
were commented because at the time there was no method defined to numerically evaluate those. This is not the case for beta()
however. Here is the evalf
method (from line 227 of ginac/inifcns_gamma.cpp
):
if (is_exactly_a<numeric>(x) && is_exactly_a<numeric>(y)) { try { return exp(lgamma(ex_to<numeric>(x))+lgamma(ex_to<numeric>(y))-lgamma(ex_to<numeric>(x+y))); } catch (const dunno &e) { } }
We'll find out when someone tries this out I suppose.
comment:6 Changed 10 years ago by
- Keywords beginner added
Changed 10 years ago by
comment:7 Changed 10 years ago by
- Status changed from new to needs_work
Added ginac wrapper for beta function. There is a problem when one argument is a float.
--Karen
sage: beta(3,2.0) ------------------------------------------------------------ Unhandled SIGSEGV: A segmentation fault occurred in Sage. This probably occurred because a *compiled* component of Sage has a bug in it (typically accessing invalid memory) or is not properly wrapped with _sig_on, _sig_off. You might want to run Sage under gdb with 'sage -gdb' to debug this. Sage will now terminate (sorry). ------------------------------------------------------------
comment:8 Changed 10 years ago by
- Keywords beginner removed
- Milestone changed from sage-5.0 to sage-4.7.1
Hi Karen,
sorry for taking so long to look at this. It seems that I forgot to check for a NULL pointer in py_float()
. attachment:trac_9130-py_float_segfault.patch fixes the segfault.
sage: from sage.functions.other import beta sage: beta(3,2.0) 0.0833333333333333
Will you have time to finish the patch?
You need to add an import statement to sage/functions/all.py
so beta
is available at the command line. The documentation also needs some care. IIRC there should be an empty line after INPUT
, OUTPUT
and EXAMPLES
. The statement
It is computed by various libraries within Sage, depending on the input type.
is too vague. We should either remove it or explain how GiNaC evaluates this (see beta_eval()
and beta_evalf()
in inifcns_gamma.cpp).
comment:9 Changed 9 years ago by
- Keywords sd35.5 added
comment:10 Changed 9 years ago by
- Cc benjaminfjones added
Changed 9 years ago by
comment:11 Changed 9 years ago by
- Status changed from needs_work to needs_review
comment:12 Changed 9 years ago by
Apply all patches in the order they were added:
- trac_9130_beta_function.patch
- trac_9130-py_float_segfault.patch
- trac_9130_beta_function.2.patch
- trac_9130_beta_function.3.patch
comment:13 Changed 9 years ago by
- Description modified (diff)
comment:14 Changed 9 years ago by
- Status changed from needs_review to needs_work
I think we discovered that the only complex inputs that the code accepts are ones where one of the parameters is equal to 1. In that case beta(1,x) = 1/x
is used to compute the result. Looks like GiNaC can't handle complex inputs at all (or perhaps complex numbers aren't being passed to GiNaC in a way it understands).
On the other hand, mpmath does support evaluation at arbitrary precision complex numbers so that could be a useful enhancement that could take place in a new ticket.
I would change the docstrings to clearly indicate that
- only real inputs are accepted (for now)
- beta(1,x) = beta(x,1) = 1/x simplification is automatically applied
comment:15 Changed 9 years ago by
- Reviewers set to Benjamin Jones
comment:16 Changed 9 years ago by
- Status changed from needs_work to needs_review
comment:17 Changed 9 years ago by
- Description modified (diff)
comment:18 Changed 9 years ago by
I uploaded a combined *single* patch that folds together the previously uploaded 5 patches applied in order. The patch looks great to me. I have one question about what happens when GinacFunction
's __call__
returns a TypeError whose message doesn't contain "cannot coerce". What gets returned in that case? (the comments in the gamma
implementation that the try / except block was borrowed from say that TypeErrors? are raised when a fast float is passed. These are then ignored in the current implementation? I haven't tested this.
I'm running make ptestlong
on Sage-4.8.alpha6 with the combined patch applied overnight; will update in the morning.
Changed 9 years ago by
comment:19 Changed 9 years ago by
- Description modified (diff)
- Milestone changed from sage-4.8 to sage-5.0
- Status changed from needs_review to needs_work
I uploaded a new patch attachment:trac_9130-py_float_segfault.take2.patch which checks for TypeError
instead of ValueError
in sage.symbolic.pynac.py_float()
. Now we can evaluate beta()
on complex values, so the __call__()
method in Karen's patches can be removed.
Changed 9 years ago by
new combined file replacement by trac_9130-py_float_segfault.take2.patch and removal of __call__()
method
comment:20 Changed 9 years ago by
- Description modified (diff)
- Status changed from needs_work to needs_review
New patch attachment:trac_9130_combined.2.patch replaces burcin's old patch with the newer one and also removes the __call__()
method in beta.
Apply only this one.
comment:21 Changed 9 years ago by
I added #12303 as a separate ticket to leave beta symbolic for exact complex inputs such as beta(2,1+5*I)
comment:22 Changed 9 years ago by
- Reviewers changed from Benjamin Jones to Benjamin Jones, Burcin Erocal
- Status changed from needs_review to needs_work
The patch looks great, I have a couple of minor comments. We can switch this to positive review when these are fixed:
- This paragraph in the docstring is not needed any more:
GiNaC is used to compute `B(p,q)`. However, complex inputs are not yet handled in general. When GiNaC raises an error on such inputs, we raise a NotImplementedError.
Perhaps the first sentence can be merged with the paragraph following that.
- The commit message for the combined patch is not very helpful. I'd like my patch to stay separate if possible.
comment:23 Changed 9 years ago by
I forgot to update -- with attachment:trac_9130_combined.patch applied to 4.8.alpha6, make ptestlong
succeeds with no failures! I'll test the updated patch now.
It should be easy to edit the newest combined patch file to remove the changes from Burcin's patches and just have those applied separately, right?
comment:24 Changed 9 years ago by
- Description modified (diff)
- Status changed from needs_work to positive_review
I uploaded a new patch which combines Karen's previous patches. This is ready to go now.
For numerical evaluation, mpmath has beta and also the generalized incomplete beta function for complex arguments. But it's probably easy to do the complete beta function directly with Ginac.
Simplification for rational arguments (beta(0.5,0.5) = pi) would be nice.
Unless someone else wants to work on this, I might have a stab at it within a couple of days.