Opened 9 years ago

Closed 8 years ago

Last modified 3 years ago

#9130 closed enhancement (fixed)

Access to beta function

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: sage-5.0.beta6
Authors: Karen Kohl, Burcin Erocal, Karl-Dieter Crisman Reviewers: Benjamin Jones, Burcin Erocal, Karl-Dieter Crisman
Report Upstream: N/A Work issues:
Branch: Commit:
Dependencies: #4498, #12507 Stopgaps:

Description (last modified by jdemeyer)

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

Attachments (15)

trac_9130_beta_function.patch (1.5 KB) - added by ktkohl 9 years ago.
trac_9130-py_float_segfault.patch (2.4 KB) - added by burcin 8 years ago.
fix segfault in py_float
trac_9130_beta_function.2.patch (4.2 KB) - added by ktkohl 8 years ago.
trac_9130_beta_function.3.patch (1.8 KB) - added by ktkohl 8 years ago.
include beta in random_tests
trac_9130_beta_function.4.patch (3.3 KB) - added by ktkohl 8 years ago.
updated comments/examples/error message text.
trac_9130_combined.patch (8.6 KB) - added by benjaminfjones 8 years ago.
combined reviewer patch folding up previously uploaded 5 patches
trac_9130-py_float_segfault.take2.patch (2.4 KB) - added by burcin 8 years ago.
trac_9130_combined.2.patch (7.4 KB) - added by ktkohl 8 years ago.
new combined file replacement by trac_9130-py_float_segfault.take2.patch and removal of __call__() method
trac_9130-beta_function.patch (5.0 KB) - added by burcin 8 years ago.
combined Karen's patches
trac_9130-reviewer.2.patch (5.7 KB) - added by kcrisman 8 years ago.
trac_9130-reviewer.patch (5.7 KB) - added by kcrisman 8 years ago.
Apply last
trac_9130-beta_function.2.patch (3.5 KB) - added by benjaminfjones 8 years ago.
trac_9130-random-tests.patch (2.5 KB) - added by benjaminfjones 8 years ago.
fixes random tests after rebasing against #4498
trac_9130-random-tests.2.patch (2.5 KB) - added by jdemeyer 8 years ago.
trac_9130-beta_function.3.patch (3.5 KB) - added by jdemeyer 8 years ago.

Download all attachments as: .zip

Change History (60)

comment:1 Changed 9 years ago by fredrik.johansson

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.

comment:2 Changed 9 years ago by kcrisman

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: Changed 9 years ago by 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 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: Changed 9 years ago by 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 in sage/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 9 years ago by burcin

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 in sage/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 9 years ago by burcin

  • Keywords beginner added

Changed 9 years ago by ktkohl

comment:7 Changed 9 years ago by ktkohl

  • 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).
------------------------------------------------------------

Changed 8 years ago by burcin

fix segfault in py_float

comment:8 Changed 8 years ago by burcin

  • Authors set to Karen T. Kohl, Burcin Erocal
  • 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 8 years ago by ktkohl

  • Keywords sd35.5 added

comment:10 Changed 8 years ago by benjaminfjones

  • Cc benjaminfjones added

Changed 8 years ago by ktkohl

comment:11 Changed 8 years ago by ktkohl

  • Status changed from needs_work to needs_review

Changed 8 years ago by ktkohl

include beta in random_tests

comment:12 Changed 8 years ago by ktkohl

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

  • Description modified (diff)

comment:14 Changed 8 years ago by benjaminfjones

  • 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

  1. only real inputs are accepted (for now)
  2. beta(1,x) = beta(x,1) = 1/x simplification is automatically applied

comment:15 Changed 8 years ago by benjaminfjones

  • Reviewers set to Benjamin Jones

Changed 8 years ago by ktkohl

updated comments/examples/error message text.

comment:16 Changed 8 years ago by ktkohl

  • Status changed from needs_work to needs_review

comment:17 Changed 8 years ago by kcrisman

  • Description modified (diff)

Changed 8 years ago by benjaminfjones

combined reviewer patch folding up previously uploaded 5 patches

comment:18 Changed 8 years ago by benjaminfjones

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

comment:19 Changed 8 years ago by burcin

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

new combined file replacement by trac_9130-py_float_segfault.take2.patch and removal of __call__() method

comment:20 Changed 8 years ago by ktkohl

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

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

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

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?

Changed 8 years ago by burcin

combined Karen's patches

comment:24 Changed 8 years ago by burcin

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

comment:25 Changed 8 years ago by burcin

  • Keywords Cernay2012 added

comment:26 Changed 8 years ago by kcrisman

Hooray!

comment:27 Changed 8 years ago by jdemeyer

  • Status changed from positive_review to needs_work

There are problems building the documentation:

/mnt/usb1/scratch/jdemeyer/merger/sage-5.0.beta4/local/lib/python2.7/site-packages/sage/functions/other.py:docstring of sage.functions.other.psi:6: WARNING: Inline literal start-string without end-string.

comment:28 Changed 8 years ago by kcrisman

  • Reviewers changed from Benjamin Jones, Burcin Erocal to Benjamin Jones, Burcin Erocal, Karl-Dieter Crisman

This is because of the addition to the reference manual of stuff that wasn't up to snuff, though almost. There was an extra colon as well, and a few other similar things.

I've fixed these things, and actually tried adding loads(dumpg()) tests as well, but it doesn't matter, because it's not testing right. It's not the computer, since other files are passing fine. It seems to be the symbolic version of beta.

Getting rid of all stuff with x and keeping

            sage: beta(1,2.0+I)
            0.400000000000000 - 0.200000000000000*I

causes

$ ../../sage -t sage/functions/other.py 
sage -t  "devel/sage-main/sage/functions/other.py"          
The doctested process was killed by signal 11
	 [6.8 s]
 
----------------------------------------------------------------------
The following tests failed:


	sage -t  "devel/sage-main/sage/functions/other.py" # Killed/crashed
Total time for all tests: 6.9 seconds

but removing that as well yields passing tests.

This seems unexpected to me.

comment:29 Changed 8 years ago by kcrisman

Sorry, ignore all that. I didn't see Burcin's other patch because I read the comments, not the description.

comment:30 Changed 8 years ago by kcrisman

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

This requires only that someone do

./sage -docbuild reference html

and look at

devel/sage-main/doc/output/html/en/reference/other.html

as well as run tests.

Changed 8 years ago by kcrisman

Changed 8 years ago by kcrisman

Apply last

comment:31 Changed 8 years ago by kcrisman

  • Authors changed from Karen T. Kohl, Burcin Erocal to Karen T. Kohl, Burcin Erocal, Karl-Dieter Crisman

To finish:

comment:32 Changed 8 years ago by benjaminfjones

  • Status changed from needs_review to positive_review

Finished! The patches apply cleanly to 5.0.beta3, I checked the docs and they look good. All tests pass (I tested sage/functions/other.py and also did make ptest for good measure).

comment:33 Changed 8 years ago by jdemeyer

  • Status changed from positive_review to needs_work

This patch conflicts with #4498. Either this one or #4498 should be rebased.

comment:34 Changed 8 years ago by benjaminfjones

  • Dependencies set to 4498
  • Description modified (diff)

comment:35 Changed 8 years ago by benjaminfjones

  • Dependencies changed from 4498 to #4498

Changed 8 years ago by benjaminfjones

Changed 8 years ago by benjaminfjones

fixes random tests after rebasing against #4498

comment:36 Changed 8 years ago by benjaminfjones

  • Description modified (diff)

comment:37 follow-ups: Changed 8 years ago by benjaminfjones

  • Status changed from needs_work to needs_review

I rebased the patches against #4498. This meant moving the changes to sage/symbolic/random_tests.py out of trac_9130_beta_function.patch and into a new patch trac_9130-random-tests.patch that reflects the differences in the random tests after #4498 was applied and beta was added.

The tests in random_test.py are getting annoying. Now #9130 must depend on #4498 and in turn my #11888 must depend on #9130, etc. To finish #11143, I need to rebase the changes to random_tests.py to all 3 of these other symbolics patches.

Would it make sense to open a new ticket for changes needed to random_tests.py after all the new symbolic functions are added that will go into sage-5.0? This new ticket would depend on all the other tickets adding new functions. This way we wouldn't have to fix the random tests like this in every ticket that adds a new function, just *once* per release.

comment:38 in reply to: ↑ 37 Changed 8 years ago by jdemeyer

Replying to benjaminfjones:

Would it make sense to open a new ticket for changes needed to random_tests.py after all the new symbolic functions are added that will go into sage-5.0? This new ticket would depend on all the other tickets adding new functions. This way we wouldn't have to fix the random tests like this in every ticket that adds a new function, just *once* per release.

I think that would be making it more complicated for yourself and for reviewers. It also means that, if just one of those tickets has a problem, none of them can be merged. But I'm not against your suggestion.

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

Replying to benjaminfjones:

I rebased the patches against #4498. This meant moving the changes to sage/symbolic/random_tests.py out of trac_9130_beta_function.patch and into a new patch trac_9130-random-tests.patch that reflects the differences in the random tests after #4498 was applied and beta was added.

Thank you for taking care of this.

The tests in random_test.py are getting annoying. Now #9130 must depend on #4498 and in turn my #11888 must depend on #9130, etc. To finish #11143, I need to rebase the changes to random_tests.py to all 3 of these other symbolics patches.

Would it make sense to open a new ticket for changes needed to random_tests.py after all the new symbolic functions are added that will go into sage-5.0? This new ticket would depend on all the other tickets adding new functions. This way we wouldn't have to fix the random tests like this in every ticket that adds a new function, just *once* per release.

At this stage I'd be happy to mark that test in random_tests.py with a #random tag so the doctesting framework ignores the output. If you open a new ticket for this and provide a patch, I promise a quick positive review. :)

comment:40 in reply to: ↑ 39 Changed 8 years ago by benjaminfjones

Replying to burcin: Thanks, marking the three "offending" doctests with a random tag is #12507.

comment:41 Changed 8 years ago by kcrisman

  • Status changed from needs_review to positive_review

Just checked one last time - yes, everything is fine! A very minor quibble is that the (why here?) added prime pi file could use a few double backticks and one or two other things, but that is immaterial in this saga.

Changed 8 years ago by jdemeyer

Changed 8 years ago by jdemeyer

comment:42 Changed 8 years ago by jdemeyer

  • Dependencies changed from #4498 to #4498, #12507
  • Description modified (diff)

Rebased again.

comment:43 Changed 8 years ago by kcrisman

I think if you would have merged this first, then #12507, there wouldn't have needed to be a rebase. Maybe the authors of #12507 should have put this as a dependency? I think we were all lax with it because it # randomized a doctest...

comment:44 Changed 8 years ago by jdemeyer

  • Merged in set to sage-5.0.beta6
  • Resolution set to fixed
  • Status changed from positive_review to closed

comment:45 Changed 3 years ago by chapoton

  • Authors changed from Karen T. Kohl, Burcin Erocal, Karl-Dieter Crisman to Karen Kohl, Burcin Erocal, Karl-Dieter Crisman
Note: See TracTickets for help on using tickets.