Opened 5 years ago

Closed 3 years ago

#16697 closed defect (fixed)

implement symbolic lower incomplete gamma function

Reported by: rws Owned by:
Priority: major Milestone: sage-7.2
Component: symbolics Keywords: gamma, incomplete, special, functions
Cc: kcrisman, nbruin Merged in:
Authors: Ralf Stephan Reviewers: Buck Evan, Paul Masson
Report Upstream: N/A Work issues:
Branch: 157b268 (Commits) Commit: 157b268d4ac45de049089e633809f89d2c3ab84e
Dependencies: #20185 Stopgaps:

Description (last modified by rws)

This is actually a defect because we leave a result from Maxima undefined:

sage: hypergeometric([1],[b],x).simplify_hypergeometric()
(b - 1)*x^(-b + 1)*e^x*gamma_greek(b - 1, x)
sage: gamma_greek
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-18-6e90901bc5cb> in <module>()
----> 1 gamma_greek

NameError: name 'gamma_greek' is not defined

https://en.wikipedia.org/wiki/Incomplete_gamma_function#Lower_incomplete_Gamma_function

Mathematica seems to have Gamma[a,z] for upper and Gamma[a,0,z] for lower; Maple seems to have upper Gamma. gamma_inc (the upper one in Sage) gets immediately converted to gamma(a,x). The symbolic functions gamma_inc==incomplete_gamma are converted and never returned to the user as expression:

sage: gamma_inc(x,x,hold=True)
gamma(x, x)
sage: incomplete_gamma(x,x,hold=True)
gamma(x, x)
sage: assume(x>0)
sage: integral(t^(s-1)*e^(-t),t,0,x)
-gamma(s, x) + gamma(s)

This ticket should deprecate "incomplete_gamma" and add the symbolic function gamma_inc_lower, leaving open the question of the global alias for and the displayed name of Function_gamma_inc.

Previous part of description:


  1. Provide all three "user input interfaces" gamma_inc, incomplete_gamma and lower_incomplete_gamma, and convert the Maxima gamma_greek to -gamma(a, x) + gamma(a)
  2. Provide all three "user input interfaces" gamma_inc, incomplete_gamma and lower_incomplete_gamma, and convert to gamma(a,x) and gamma(a,0,x)n on output
  3. Provide all three "user input interfaces" gamma_inc, incomplete_gamma and lower_incomplete_gamma, and have incomplete_gamma and lower_incomplete_gamma as result instead of `gamma(...)'
  4. Change the user interface completely to lower_incomplete_gamma and upper_incomplete_gamma
  5. Change the user interface completely to gamma_inc_lower and gamma_inc_upper
  6. Change the user interface completely to gamma(a,x) and gamma(a,0,x)

Attachments (2)

16697-gammainc-alt-patch (2.0 KB) - added by rws 5 years ago.
see comment
ptest.log (652.5 KB) - added by buck 4 years ago.
ptest results

Download all attachments as: .zip

Change History (82)

comment:1 Changed 5 years ago by rws

  • Cc kcrisman nbruin added

comment:2 Changed 5 years ago by nbruin

We don't need all of gamma, gamma_inc, incomplete_gamma as symbolic function if they all mean the same. Some of them are probably there for user convenience or because of legacy. I think we want only one spelling for gamma_upper, whatever it may be. The gamma(a,x) convention seems fairly prominent in computer algebra, but is a poor translation of math notation, where the difference between upper and lower is made with case (and indeed, the complete gamma function tends to be capital gamma. I don't think I've ever seen different).

Because of tab completion, you definitely want a binding available that starts with gamma. I'd be fine with gamma_lower but if you want to go with gamma_lower_incomplete, that's fine with me too. So, I'd propose

  1. Maintain status quo with respect to upper incomplete gamma, i.e., don't change sage's behaviour for gamma(a,x) and gamma_inc(a,x). Just implement a symbolic function gamma_lower (or modified spelling) to correspond to maxima's gamma_greek.
Last edited 5 years ago by nbruin (previous) (diff)

comment:3 Changed 5 years ago by rws

I don't underswtand. With "leave as is" you would keep gamma(a,x) as returned expression for the upper function? Or would you make gamma_inc fully symbolic in that it gets no longer converted to gamma(a,x)?

comment:4 Changed 5 years ago by rws

The reason why gamma_inc(x,y) gets output as gamma(x,y) is not conversion as I thought but simple naming super in the constructor. It would be just as easy to make that gamma_inc and change a few doctests such that then:

         sage: f = exp_integral_e(-1,x)
         sage: f.integrate(x)
-        Ei(-x) - gamma(-1, x)
+        Ei(-x) - gamma_inc(-1, x)

             sage: gamma_inc(3,2)
-            gamma(3, 2)
+            gamma_inc(3, 2)

             sage: gamma(3,2)
             gamma_inc(3, 2)

             sage: gamma(x,0)
             gamma(x)

             sage: gamma(x, 0, hold=True)
-            gamma(x, 0)
+            gamma_inc(x, 0)

I'm appending the patch that enables this here, so it is preserved.

Changed 5 years ago by rws

see comment

comment:5 Changed 5 years ago by rws

  • Branch set to u/rws/implement_symbolic_lower_incomplete_gamma_function

comment:6 Changed 5 years ago by git

  • Commit set to 93c2f94cb7e1ca7651a8a8de907bad09ae82428e

Branch pushed to git repo; I updated commit sha1. New commits:

93c2f9416697: deprecate "incomplete_gamma"

comment:7 Changed 5 years ago by rws

  • Description modified (diff)

comment:8 Changed 5 years ago by git

  • Commit changed from 93c2f94cb7e1ca7651a8a8de907bad09ae82428e to 5c1dd67bfa51d7c871b055c3272a94af9fdb68f2

Branch pushed to git repo; I updated commit sha1. New commits:

5c1dd6716697: implement gamma_inc_lower

comment:9 Changed 5 years ago by rws

  • Status changed from new to needs_review

comment:10 Changed 5 years ago by git

  • Commit changed from 5c1dd67bfa51d7c871b055c3272a94af9fdb68f2 to d8d7ee29d7d0f8cf9f884c4403fe5c36cdf3336d

Branch pushed to git repo; I updated commit sha1. New commits:

d8d7ee216697: use mpmath as default evalf algorithm

comment:11 follow-up: Changed 5 years ago by nbruin

On sage-support Rafael Greenblatt reports

sage: (incomplete_gamma(x,y).diff(x)).simplify()
TypeError: ECL says: Error executing code in Maxima: gamma: wrong number of arguments.

The problem probably arises because the function ends up separated from its arguments:

sage: incomplete_gamma(x,y).diff(x)
D[0](gamma)(x, y)

which might confuse the translator (the behaviour is consistent with just going by strings). The changes here might improve the situation. It's worth checking.

The following of course works because it does NOT primarily look at strings:

sage: from sage.interfaces.maxima_lib import *
sage: maxima_calculus(sr_to_max(incomplete_gamma(x,y).diff(x)))
'diff(gamma_incomplete(x,y),x,1)

comment:12 in reply to: ↑ 11 ; follow-up: Changed 5 years ago by rws

Replying to nbruin:

On sage-support Rafael Greenblatt reports

sage: (incomplete_gamma(x,y).diff(x)).simplify()
TypeError: ECL says: Error executing code in Maxima: gamma: wrong number of arguments.

The problem probably arises because the function ends up separated from its arguments:

sage: incomplete_gamma(x,y).diff(x)
D[0](gamma)(x, y)

which might confuse the translator (the behaviour is consistent with just going by strings). The changes here might improve the situation. It's worth checking.

Does the same with the patch. Maybe include the fix here? Is it simple?

comment:13 in reply to: ↑ 12 Changed 5 years ago by nbruin

Replying to rws:

Does the same with the patch. Maybe include the fix here? Is it simple?

Might not be too bad. I think the problem is in sage.symbolic.expression_conversions.InterfaceInit.derivative (line 515 or so). You can see there that the string representation of a derivative expression gets assembled by using f.name(). Replacing that by f._maxima_init_() might already do a better job:

sage: gamma(x,y).operator().name()
'gamma'
sage: gamma(x).operator().name()
'gamma'
sage: gamma(x,y).operator()._maxima_init_()
'gamma_incomplete'
sage: gamma(x).operator()._maxima_init_()
'gamma'

comment:14 follow-up: Changed 5 years ago by rws

That gives:

sage: (incomplete_gamma(x,y).diff(x)).simplify()
D[0](gamma)(x, y)

comment:15 in reply to: ↑ 14 Changed 5 years ago by nbruin

Replying to rws:

That gives:

sage: (incomplete_gamma(x,y).diff(x)).simplify()
D[0](gamma)(x, y)

Cool! so it IS easy. Be sure to also test (I haven't checked the answer is correct, but at least it shouldn't raise an error):

sage: (incomplete_gamma(x,x+1).diff(x)).simplify().simplify()
-(x + 1)^(x - 1)*e^(-x - 1) + D[0](gamma)(x, x + 1)

to check that both occurrences of f.name() have been replaced. I check the rest of that file and didn't find any further suspicious uses of name so after this, our maxima interface is probably perfect :-).

comment:16 Changed 5 years ago by nbruin

See #16785 for a more complete fix.

comment:17 Changed 5 years ago by git

  • Commit changed from d8d7ee29d7d0f8cf9f884c4403fe5c36cdf3336d to a73c5e374b1ee694ba20af22654767cdf441c126

Branch pushed to git repo; I updated commit sha1. New commits:

862798bMerge branch 'develop' into t/16697/implement_symbolic_lower_incomplete_gamma_function
705c34e16697: fix a related bug
a73c5e316697: fix previous patch and a doctest

comment:18 Changed 5 years ago by rws

I also had to fix a doctest in integration where the given result didn't have enough precision because it came from maxima. Now the default is mpmath and that uncovered it.

comment:19 Changed 5 years ago by git

  • Commit changed from a73c5e374b1ee694ba20af22654767cdf441c126 to 1fc19e73cb83d39a17d2055ec9b661f1fb401dfd

Branch pushed to git repo; I updated commit sha1. New commits:

7558a7c16697: revert fix in favor of merge of 16785
9c0f66dtrac #16785: improve differential operator translation to maxima
1fc19e7Merge branch 't/16785/ticket/16785' into t/16697/implement_symbolic_lower_incomplete_gamma_function

comment:20 Changed 5 years ago by rws

  • Merged in set to #16785

comment:21 Changed 5 years ago by vbraun_spam

  • Milestone changed from sage-6.3 to sage-6.4

comment:22 Changed 5 years ago by kcrisman

  • Merged in #16785 deleted

Sorry, we should make this clearer - that field is one currently not used but which should be (!!!) that was for saying which version/beta of Sage a given ticket was merged in. I agree that might not be clear given the word "merge" meaning something else in git!

comment:23 Changed 5 years ago by rws

So, is there anything missing?

comment:24 Changed 5 years ago by rws

  • Authors set to Ralf Stephan

comment:25 Changed 5 years ago by git

  • Commit changed from 1fc19e73cb83d39a17d2055ec9b661f1fb401dfd to 9033be83a52eaabe0c0f896193923a040837e0cf

Branch pushed to git repo; I updated commit sha1. New commits:

9033be8Merge branch 'develop' into t/16697/implement_symbolic_lower_incomplete_gamma_function

comment:26 Changed 5 years ago by git

  • Commit changed from 9033be83a52eaabe0c0f896193923a040837e0cf to 05195e75302bb598d61d724d8635b08cc1712846

Branch pushed to git repo; I updated commit sha1. New commits:

05195e7Merge branch 'develop' into t/16697/implement_symbolic_lower_incomplete_gamma_function

comment:27 Changed 5 years ago by jdemeyer

Just to let you know: this conflicts with #17130.

comment:28 Changed 5 years ago by kcrisman

See also the information here, which seems to be more information than the documentation of Maxima gives.


Implementation of the Incomplete Gamma function

 New Maxima User function: gamma_incomplete(a,z)

 The following features are implemented:

 - Evaluation for real and complex numbers in double float and
   bigfloat precision
 - Special values for gamma_incomplete(a,0) and gamma_incomplete(a,inf)
 - When $gamma_expand T expand the following expressions:
   gamma_incomplete(0,z)
   gamma_incomplete(n+1/2)
   gamma_incomplete(1/2-n)
   gamma_incomplete(n,z)
   gamma_incomplete(-n,z)
   gamma_incomplete(a+n,z)
   gamma_incomplete(a-n,z)
 - Mirror symmetry
 - Derivative wrt the arguments a and z

--------------------------------------------------------------------------------
Implementation of the Generalized Incomplete Gamma function

 New Maxima User function: gamma_incomplete_generalized(a,z1,z2)

 The following features are implemented:

 - Evaluation for real and complex numbers in double float and
   bigfloat precision
 - Special values for:
   gamma_incomplete_generalized(a,z1,0)
   gamma_incomplete_generalized(a,0,z2),
   gamma_incomplete_generalized(a,z1,inf)
   gamma_incomplete_generalized(a,inf,z2)
   gamma_incomplete_generalized(a,0,inf)
   gamma_incomplete_generalized(a,x,x)
 - When $gamma_expand T and n an integer expand
   gamma_incomplete_generalized(a+n,z1,z2)
 - Implementation of Mirror symmetry
 - Derivative wrt the arguments a, z1 and z2

--------------------------------------------------------------------------------
Implementation of the Regularized Incomplete Gamma function

 New Maxima User function: gamma_incomplete_regularized(a,z)

 The following features are implemented:

 - Evaluation for real and complex numbers in double float and
   bigfloat precision
 - Special values for:
   gamma_incomplete_regularized(a,0)
   gamma_incomplete_regularized(0,z)
   gamma_incomplete_regularized(a,inf)
 - When $gamma_expand T and n a positive integer expansions for
   gamma_incomplete_regularized(n+1/2,z)
   gamma_incomplete_regularized(1/2-n,z)
   gamma_incomplete_regularized(n,z)
   gamma_incomplete_regularized(a+n,z)
   gamma_incomplete_regularized(a-n,z)
 - Derivative wrt the arguments a and z
 - Implementation of Mirror symmetry

comment:29 Changed 5 years ago by rws

We don't have gamma_incomplete_regularized either. This would be a different ticket.

comment:30 follow-up: Changed 5 years ago by rws

I'm going back to Pari as default since mpmath has similar problems as earlier Pari, see #17328 and https://github.com/fredrik-johansson/mpmath/issues/301

comment:31 Changed 5 years ago by git

  • Commit changed from 05195e75302bb598d61d724d8635b08cc1712846 to 6f59f3a6133b9654f349a1228a6bc9fe9c90c944

Branch pushed to git repo; I updated commit sha1. New commits:

6d10729Simplify numerical evaluation of BuiltinFunctions
b6e1ed4Merge remote-tracking branches 'origin/u/jdemeyer/ticket/17131' and 'origin/u/jdemeyer/ticket/17133' into ticket/17130
382f97aMerge branch 'u/jdemeyer/ticket/17130' of trac.sagemath.org:sage into 6.5beta1
726598917130: reviewer's patch: fix typo
c47dbd4Fix Trac #17328 again in a better way
a486db2Call the factorial() method of an Integer
9d3cbbdFix numerical noise
abab222Fix more numerical noise
a383a5bMerge branch 'u/jdemeyer/ticket/17130' of trac.sagemath.org:sage into t/16697/implement_symbolic_lower_incomplete_gamma_function
6f59f3a16697: adaptations due to trac 17130

comment:32 Changed 5 years ago by rws

  • Dependencies set to #17130

comment:33 Changed 5 years ago by chapoton

  • Branch changed from u/rws/implement_symbolic_lower_incomplete_gamma_function to public/ticket/16697
  • Commit changed from 6f59f3a6133b9654f349a1228a6bc9fe9c90c944 to 202b202c0fd19a4848997e7c18895a216d9ad57b

rebased on latest 6.5.beta3

also added one missing trac role


New commits:

202b202Merge branch 'u/rws/implement_symbolic_lower_incomplete_gamma_function' of trac.sagemath.org:sage into 6.5.b3

comment:34 Changed 5 years ago by git

  • Commit changed from 202b202c0fd19a4848997e7c18895a216d9ad57b to f8fcbc25f545e412f6c1cf18d58e7f50f4e44cb3

Branch pushed to git repo; I updated commit sha1. New commits:

f8fcbc2trac #16697 one more trac role missing was added

comment:35 follow-up: Changed 5 years ago by kcrisman

Frédéric, any comments on review, or was this just a drive-by improvement?

Does the following (from the description) now come out better?

sage: hypergeometric([1],[b],x).simplify_hypergeometric()

Could we add the gamma(a,0,x) interface for the lower gamma? Should we?

comment:36 in reply to: ↑ 35 ; follow-up: Changed 5 years ago by rws

Replying to kcrisman:

Does the following (from the description) now come out better?

sage: hypergeometric([1],[b],x).simplify_hypergeometric()

(b - 1)*x^(-b + 1)*e^x*gamma_inc_lower(b - 1, x)

Could we add the gamma(a,0,x) interface for the lower gamma? Should we?

A separate ticket, please?

comment:37 in reply to: ↑ 36 ; follow-up: Changed 5 years ago by kcrisman

sage: hypergeometric([1],[b],x).simplify_hypergeometric()

(b - 1)*x^(-b + 1)*e^x*gamma_inc_lower(b - 1, x)

Nice.

Could we add the gamma(a,0,x) interface for the lower gamma? Should we?

A separate ticket, please?

Seems quite reasonable. Do you agree this would be a useful addition? I only mention it because it seems that at least a couple 'competitors' support it, so presumably it is somewhat of a standard and not a surprise if Sage implements it.

comment:38 in reply to: ↑ 37 Changed 5 years ago by rws

Replying to kcrisman:

Could we add the gamma(a,0,x) interface for the lower gamma? Should we?

A separate ticket, please?

Seems quite reasonable. Do you agree this would be a useful addition? I only mention it because it seems that at least a couple 'competitors' support it, so presumably it is somewhat of a standard and not a surprise if Sage implements it.

Any UI change that reduces surprise is fine with me. (Actually, I have a GUI development history and have read and care much about this subject)

comment:39 Changed 5 years ago by kcrisman

For the record, I know nearly nothing about UI development! It was only a question. #17722

I'm going to hold off on looking at this in more detail in case chapoton already did but a brief glance looks good.

comment:40 Changed 5 years ago by kcrisman

Note sympy apparently implements this as lowergamma (?).

comment:41 Changed 5 years ago by rws

  • Milestone changed from sage-6.4 to sage-6.6
  • Status changed from needs_review to needs_work

Thanks, I was completely unaware. I think I will use their eval. The diff cannot be used in completeness because we don't know about the Meijer G function (EDIT: #17970), but I can steal half of it.

Last edited 5 years ago by rws (previous) (diff)

comment:42 Changed 5 years ago by git

  • Commit changed from f8fcbc25f545e412f6c1cf18d58e7f50f4e44cb3 to 5f4a7fbd8598ef3fb22d8cb61d199dfffea23a8c

Branch pushed to git repo; I updated commit sha1. New commits:

9cc532cMerge branch 'develop' into t/16697/public/ticket/16697
5f4a7fb16697: eval using sympy; diff; doctests

comment:43 Changed 5 years ago by rws

  • Dependencies changed from #17130 to #17130, sympy-0.7.7
  • Milestone changed from sage-6.6 to sage-pending
  • Status changed from needs_work to needs_review

Three doctests fail now because we depend on sympy.erf._sage_() which will be in Sympy-0.7.7. So this is pending.

comment:44 follow-up: Changed 5 years ago by buck

Is there a ticket for sympy-0.7.7?

comment:45 in reply to: ↑ 44 Changed 5 years ago by rws

Replying to buck:

Is there a ticket for sympy-0.7.7?

There will be as soon it is released!

comment:46 Changed 5 years ago by rws

The diff link of this ticket shows not the correct diff but if you git trac checkout 16697 you get the right branch.

comment:47 follow-up: Changed 5 years ago by buck

I tried to review this (per ticket:18210#comment:18) but quickly got confused. (This would be a lot easier on github.)

  1. Could one of you point me to the doc pertaining to _sage_? I did look, but the special sage functions page doesn't mention it.
  2. There are some doctests regarding precision that have been obliterated. Was this intentional? Why?
  3. It seems like this new function (gamma_inc_lower) would be deprecated as soon as #17722 is implemented. Don't we want to take a minute to figure out a design that works for both issues?

I think if this patch implemented lower-gamma as a special case of a single general gamma function we would end up with less if statements / edge cases / opportunities for bugs.

Edit: ... end up with *less* bugs

Last edited 5 years ago by buck (previous) (diff)

comment:48 in reply to: ↑ 47 Changed 5 years ago by rws

Replying to buck:

  1. Could one of you point me to the doc pertaining to _sage_?

https://github.com/sagemath/sage/blob/master/src/sage/structure/sage_object.pyx#L557 This is overridden where necessary so please see the resp. class. Probably you want to know about interfaces: https://github.com/sagemath/sage/blob/master/src/sage/interfaces/interface.py#L803 http://sagemath.org/doc/reference/interfaces/index.html

  1. There are some doctests regarding precision that have been obliterated. Was this intentional? Why?

Oh that's true, I don't know why that got ditched.

  1. It seems like this new function (gamma_inc_lower) would be deprecated as soon as #17722 is implemented. Don't we want to take a minute to figure out a design that works for both issues?

There should be no problem to support both interfaces, and without code duplication. We already have

sage: incomplete_gamma(1,2)
e^(-2)
sage: gamma(1,2)
e^(-2)
sage: gamma(2)
1

I think if this patch implemented lower-gamma as a special case of a single general gamma function we would end up with less if statements / edge cases / opportunities for bugs.

Well that's the task for #17722, right. To have many doctests will help a good deal.

Note that even if you review this ticket (please add your real name to the ticket head in case) we would need to wait for Sympy-0.7.7 to be released. But a review of the code and the math would leave only a make ptestlong to finish it all.

EDIT: typos

Last edited 5 years ago by rws (previous) (diff)

comment:49 Changed 5 years ago by buck

  • Reviewers set to Buck Evan
  • Status changed from needs_review to needs_work

comment:50 Changed 4 years ago by buck

I left this as "needs work" because of your "Oh that's true, I don't know why that got ditched."

Could you merge this up with master and I'll attempt to review again?

Also, please show me how to run all the relevant tests, please?

I'm surprised sympy hasn't cut a release yet, still. If I'm counting correctly, it's been eight months since their last release.

comment:51 Changed 4 years ago by buck

Merged and attempted to run the doctests, and got this result:

./sage -t src/sage/functions/other.py                                                                    ⬆ ✱ ◼
too few successful tests, not using stored timings
Running doctests with ID 2015-07-26-17-39-51-10ac5315.
Git branch: master
Using --optional=mpir,python2,sage,scons
Doctesting 1 file.
sage -t src/sage/functions/other.py
**********************************************************************
File "src/sage/functions/other.py", line 1076, in sage.functions.other.Function_gamma_inc_lower._eval_
Failed example:
    gamma_inc_lower(1/2,2)
Exception raised:
    Traceback (most recent call last):
      File "/home/buck/trees/mine/sage/local/lib/python2.7/site-packages/sage/doctest/forker.py", line 496, in _run
        self.compile_and_execute(example, compiler, test.globs)
      File "/home/buck/trees/mine/sage/local/lib/python2.7/site-packages/sage/doctest/forker.py", line 858, in compile_and_execute
        exec(compiled, globs)
      File "<doctest sage.functions.other.Function_gamma_inc_lower._eval_[2]>", line 1, in <module>
        gamma_inc_lower(Integer(1)/Integer(2),Integer(2))
      File "sage/symbolic/function.pyx", line 994, in sage.symbolic.function.BuiltinFunction.__call__ (/home/buck/trees/mine/sage/src/build/cythonized/sage/symbolic/function.cpp:10875)
        res = super(BuiltinFunction, self).__call__(
      File "sage/symbolic/function.pyx", line 507, in sage.symbolic.function.Function.__call__ (/home/buck/trees/mine/sage/src/build/cythonized/sage/symbolic/function.cpp:6887)
        res = g_function_eval2(self._serial, (<Expression>args[0])._gobj,
      File "sage/symbolic/pynac.pyx", line 185, in sage.symbolic.pynac.pyExpression_to_ex (/home/buck/trees/mine/sage/src/build/cythonized/sage/symbolic/pynac.cpp:3934)
        t = ring.SR.coerce(res)
      File "sage/structure/parent.pyx", line 1310, in sage.structure.parent.Parent.coerce (/home/buck/trees/mine/sage/src/build/cythonized/sage/structure/parent.c:10682)
        cpdef coerce(self, x):
      File "sage/structure/parent.pyx", line 1339, in sage.structure.parent.Parent.coerce (/home/buck/trees/mine/sage/src/build/cythonized/sage/structure/parent.c:10628)
        return (<map.Map>mor)._call_(x)
      File "sage/symbolic/ring.pyx", line 927, in sage.symbolic.ring.UnderscoreSageMorphism._call_ (/home/buck/trees/mine/sage/src/build/cythonized/sage/symbolic/ring.cpp:10428)
        return self.codomain()(a._sage_())
      File "/home/buck/trees/mine/sage/local/lib/python2.7/site-packages/sympy/core/mul.py", line 1467, in _sage_
        s *= x._sage_()
    AttributeError: 'erf' object has no attribute '_sage_'
**********************************************************************
File "src/sage/functions/other.py", line 1078, in sage.functions.other.Function_gamma_inc_lower._eval_
Failed example:
    gamma_inc_lower(1/2,1)
Exception raised:
    Traceback (most recent call last):
      File "/home/buck/trees/mine/sage/local/lib/python2.7/site-packages/sage/doctest/forker.py", line 496, in _run
        self.compile_and_execute(example, compiler, test.globs)
      File "/home/buck/trees/mine/sage/local/lib/python2.7/site-packages/sage/doctest/forker.py", line 858, in compile_and_execute
        exec(compiled, globs)
      File "<doctest sage.functions.other.Function_gamma_inc_lower._eval_[3]>", line 1, in <module>
        gamma_inc_lower(Integer(1)/Integer(2),Integer(1))
      File "sage/symbolic/function.pyx", line 994, in sage.symbolic.function.BuiltinFunction.__call__ (/home/buck/trees/mine/sage/src/build/cythonized/sage/symbolic/function.cpp:10875)
        res = super(BuiltinFunction, self).__call__(
      File "sage/symbolic/function.pyx", line 507, in sage.symbolic.function.Function.__call__ (/home/buck/trees/mine/sage/src/build/cythonized/sage/symbolic/function.cpp:6887)
        res = g_function_eval2(self._serial, (<Expression>args[0])._gobj,
      File "sage/symbolic/pynac.pyx", line 185, in sage.symbolic.pynac.pyExpression_to_ex (/home/buck/trees/mine/sage/src/build/cythonized/sage/symbolic/pynac.cpp:3934)
        t = ring.SR.coerce(res)
      File "sage/structure/parent.pyx", line 1310, in sage.structure.parent.Parent.coerce (/home/buck/trees/mine/sage/src/build/cythonized/sage/structure/parent.c:10682)
        cpdef coerce(self, x):
      File "sage/structure/parent.pyx", line 1339, in sage.structure.parent.Parent.coerce (/home/buck/trees/mine/sage/src/build/cythonized/sage/structure/parent.c:10628)
        return (<map.Map>mor)._call_(x)
      File "sage/symbolic/ring.pyx", line 927, in sage.symbolic.ring.UnderscoreSageMorphism._call_ (/home/buck/trees/mine/sage/src/build/cythonized/sage/symbolic/ring.cpp:10428)
        return self.codomain()(a._sage_())
      File "/home/buck/trees/mine/sage/local/lib/python2.7/site-packages/sympy/core/mul.py", line 1467, in _sage_
        s *= x._sage_()
    AttributeError: 'erf' object has no attribute '_sage_'
**********************************************************************
File "src/sage/functions/other.py", line 1092, in sage.functions.other.Function_gamma_inc_lower._eval_
Failed example:
    gamma_inc_lower(9/2,37/7)
Exception raised:
    Traceback (most recent call last):
      File "/home/buck/trees/mine/sage/local/lib/python2.7/site-packages/sage/doctest/forker.py", line 496, in _run
        self.compile_and_execute(example, compiler, test.globs)
      File "/home/buck/trees/mine/sage/local/lib/python2.7/site-packages/sage/doctest/forker.py", line 858, in compile_and_execute
        exec(compiled, globs)
      File "<doctest sage.functions.other.Function_gamma_inc_lower._eval_[10]>", line 1, in <module>
        gamma_inc_lower(Integer(9)/Integer(2),Integer(37)/Integer(7))
      File "sage/symbolic/function.pyx", line 994, in sage.symbolic.function.BuiltinFunction.__call__ (/home/buck/trees/mine/sage/src/build/cythonized/sage/symbolic/function.cpp:10875)
        res = super(BuiltinFunction, self).__call__(
      File "sage/symbolic/function.pyx", line 507, in sage.symbolic.function.Function.__call__ (/home/buck/trees/mine/sage/src/build/cythonized/sage/symbolic/function.cpp:6887)
        res = g_function_eval2(self._serial, (<Expression>args[0])._gobj,
      File "sage/symbolic/pynac.pyx", line 185, in sage.symbolic.pynac.pyExpression_to_ex (/home/buck/trees/mine/sage/src/build/cythonized/sage/symbolic/pynac.cpp:3934)
        t = ring.SR.coerce(res)
      File "sage/structure/parent.pyx", line 1310, in sage.structure.parent.Parent.coerce (/home/buck/trees/mine/sage/src/build/cythonized/sage/structure/parent.c:10682)
        cpdef coerce(self, x):
      File "sage/structure/parent.pyx", line 1339, in sage.structure.parent.Parent.coerce (/home/buck/trees/mine/sage/src/build/cythonized/sage/structure/parent.c:10628)
        return (<map.Map>mor)._call_(x)
      File "sage/symbolic/ring.pyx", line 927, in sage.symbolic.ring.UnderscoreSageMorphism._call_ (/home/buck/trees/mine/sage/src/build/cythonized/sage/symbolic/ring.cpp:10428)
        return self.codomain()(a._sage_())
      File "/home/buck/trees/mine/sage/local/lib/python2.7/site-packages/sympy/core/add.py", line 732, in _sage_
        s += x._sage_()
      File "/home/buck/trees/mine/sage/local/lib/python2.7/site-packages/sympy/core/mul.py", line 1467, in _sage_
        s *= x._sage_()
    AttributeError: 'erf' object has no attribute '_sage_'
**********************************************************************
1 item had failures:
   3 of  12 in sage.functions.other.Function_gamma_inc_lower._eval_
    [509 tests, 3 failures, 3.90 s]
----------------------------------------------------------------------
sage -t src/sage/functions/other.py  # 3 doctests failed
----------------------------------------------------------------------
Total time for all tests: 4.1 seconds
    cpu time: 2.5 seconds
    cumulative wall time: 3.9 seconds

comment:52 Changed 4 years ago by buck

I realized I needed to upgrade sympy for this branch to work. After that, the functions.other doctests pass, but a full ptest shows some failures that seem unrelated, but it's hard for me to be sure. I've attached the full result.

Changed 4 years ago by buck

ptest results

comment:53 Changed 4 years ago by buck

There's something not quite right with this branch. I get gamma_greek if I use combine().

sage: e
gamma(x, y) + gamma_inc_lower(x, y)
sage: e.combine()
gamma(x, y) + gamma_greek(x, y)
sage: e(x=1, y=2)
1
sage: e.combine()(x=1, y=2)
e^(-2) + gamma_greek(1, 2)

comment:54 Changed 4 years ago by rws

  • Branch changed from public/ticket/16697 to u/rws/16697-1

comment:55 Changed 4 years ago by rws

  • Commit changed from 5f4a7fbd8598ef3fb22d8cb61d199dfffea23a8c to f3be6f16908734f6c6e752578653a3ca09189011
  • Dependencies changed from #17130, sympy-0.7.7 to sympy-0.7.7
  • Status changed from needs_work to needs_review

Cannot confirm your doctest failures (except those where sympy-0.7.7 is needed), nor the combine problem. It looks like you tested with an older Sage the doctests of a newer Sage.

I have squashed all commits into one, and uploaded a new branch based on Sage-6.8. In this branch the three missing doctests were added as well.

Really, a full test should only be done with sympy-0.7.7 but if you think the rest is fine(?) then the missing long test eventually is done quickly by the patchbot.


New commits:

f3be6f116697: implement symbolic lower incomplete gamma function

comment:56 follow-up: Changed 4 years ago by buck

I did the full test with sympy-0.7.7, only after I got your new doctests in functions.other passing.

I reproduced the combine-gamma-greek problem with sage -br, which I would think should avoid any version-mismatch issues no? (Apparently the maxima writer forgot that uppercase gamma is also greek?)

That said, I can't reproduce it anymore after merging your branch and running another sage -br.

Should this not come out as gamma_inc_lower?

./sage -br
sage: maxima('gamma_greek(x, y)')
gamma_greek(x,y)

I do see that _sage_ of the same expression yields gamma_inc_lower. I suppose I don't know enough about the coersion going on here. I did try looking in the docs but didn't find anything that helped me understand this particular aspect.

sage: maxima('gamma_greek(x, y)')
gamma_greek(x,y)

sage: maxima('gamma_greek(x, y)')._sage_()
gamma_inc_lower(x, y)

sage: maxima('gamma_greek(x, y)').simplify()
simplify(gamma_greek(x,y))

sage: maxima('gamma_greek(x, y)')._sage_().simplify()
gamma_inc_lower(x, y)

comment:57 in reply to: ↑ 56 Changed 4 years ago by rws

Replying to buck:

I did the full test with sympy-0.7.7, only after I got your new doctests in functions.other passing.

How did you do that. Are there release tarballs of sympy-0.7.7? Note that having it on your machine does not mean Sage will use it. The version used is that in the upstream directory of Sage.

comment:58 follow-up: Changed 4 years ago by buck

I installed their master branch to the local/ tree. Putting local/bin on the $PATH makes it close enough to a virtualenv that pip will Just Work.

I needed this to get the doctests working as well, so I presumed it was also your procedure. If not, how were you getting your tests to pass?

comment:59 in reply to: ↑ 58 Changed 4 years ago by rws

Replying to buck:

I needed this to get the doctests working as well, so I presumed it was also your procedure. If not, how were you getting your tests to pass?

Never said I did. The ticket has pending and dependency set as well. But thanks anyway 8) So, except from a final test, would you say your review is positive?

comment:60 Changed 4 years ago by buck

There's doctest coverage missing (for the deprecated function), and this patch doesn't fix the numerical instability in incomplete gamma as promised elsewhere.

sage: gamma_inc_lower(25, 14.5).n()
-6.63538851954338e22
sage: gamma_inc_lower(25, 14.5).n(algorithm='mpmath')
-6.63538851954338e22
sage: gamma_inc_lower(25, 14.5).n(algorithm='pari')
-6.63538851954338e22

comment:61 Changed 4 years ago by buck

Of note, gamma (even incomplete) is strictly positive in the reals; it's an integral of a strictly positive function.

https://en.wikipedia.org/wiki/Incomplete_gamma_function#Definition

comment:62 Changed 4 years ago by buck

Wolfram gives the value of the same as 4.7e21, which seems correct.

http://www.wolframalpha.com/input/?i=Gamma%5B25%2C+0%2C+14.5%5D

comment:63 follow-up: Changed 4 years ago by buck

Actually, if I ask pari directly, it gets the correct value, so maybe the algorithm argument is broken?

sage: pari('gamma(25) - incgam(25, 14.5)')
4.71197173824555 E21
sage: pari('incgamc(25, 14.5)')
4.71197173824555 E21

comment:64 follow-up: Changed 4 years ago by buck

Same story if I ask mpmath directly. I can't imagine what's wrong.

sage: call(mpmath.gammainc, 25, 0, 14.5)
4.71197173824555e21

comment:65 in reply to: ↑ 64 Changed 4 years ago by rws

Replying to buck:

Same story if I ask mpmath directly. I can't imagine what's wrong.

'algorithm' isn't propagated to evalf so we always get the buggy pari values. As we would have to switch the default anyway as soon as a later pari version fixes the bug, I'll now simply disable pari. We'll have to live with mpmath's 53 bits for now.

comment:66 in reply to: ↑ 63 Changed 4 years ago by rws

Replying to buck:

Actually, if I ask pari directly, it gets the correct value, so maybe the algorithm argument is broken?

sage: pari('gamma(25) - incgam(25, 14.5)')
4.71197173824555 E21
sage: pari('incgamc(25, 14.5)')
4.71197173824555 E21

This works because it uses the gp interface (the libpari might be broken on our side). Ah, that's two different bugs. However, even if it is in our interface to libpari and we fix this, or if we use the gp interface, there is still the unfixed pari bug at http://pari.math.u-bordeaux.fr/cgi-bin/bugreport.cgi?bug=1689PARI/GP

comment:67 Changed 4 years ago by git

  • Commit changed from f3be6f16908734f6c6e752578653a3ca09189011 to 86618379012a463b10f0d392f57403755300236a

Branch pushed to git repo; I updated commit sha1. New commits:

866183716697: change algorithm; add, fix doctests

comment:68 Changed 4 years ago by rws

Actually the 53-bit restriction applies only to larger values, so all is not lost.

comment:69 in reply to: ↑ 30 Changed 4 years ago by rws

Replying to rws:

I'm going back to Pari as default since mpmath has similar problems as earlier Pari, see #17328 and https://github.com/fredrik-johansson/mpmath/issues/301

The mpmath issue is now fixed as well so we have both options, mpmath or Pari/GP.

comment:70 Changed 4 years ago by fredrik.johansson

You could also compute the incomplete gamma function with Arb.

comment:71 Changed 4 years ago by rws

  • Branch changed from u/rws/16697-1 to u/rws/16697-2

comment:72 Changed 4 years ago by rws

  • Commit changed from 86618379012a463b10f0d392f57403755300236a to c6be42fd29966c4b824010c998820c37169c89c0
  • Dependencies changed from sympy-0.7.7 to #20185
  • Milestone changed from sage-pending to sage-7.2

New commits:

c6be42fMerge branch 'u/rws/16697-1' of git://trac.sagemath.org/sage into i16697

comment:73 Changed 4 years ago by chapoton

  • Status changed from needs_review to needs_work

does not apply

comment:74 Changed 4 years ago by rws

Wanna review?

comment:75 Changed 3 years ago by git

  • Commit changed from c6be42fd29966c4b824010c998820c37169c89c0 to 157b268d4ac45de049089e633809f89d2c3ab84e

Branch pushed to git repo; I updated commit sha1. New commits:

1f1c367Merge branch 'develop' into t/16697/16697-2
157b268remove some proposed but obsolete doctests

comment:76 Changed 3 years ago by rws

  • Status changed from needs_work to needs_review

comment:77 Changed 3 years ago by paulmasson

Current 7.3.beta5 Sage builds fine with these modifications. Doctests all pass. Documentation builds.

Random numeric tests accurate. Function plots. One curious difference from gamma_inc:

sage: gamma_inc_lower(a,x).diff(x)
x^(a - 1)*e^(-x)
sage: gamma_inc(a,x).diff(x)
D[1](gamma)(a, x)

Looks good to me. Sympy already updated: any other reason not to merge?

comment:78 Changed 3 years ago by paulmasson

  • Reviewers changed from Buck Evan to Buck Evan, Paul Masson
  • Status changed from needs_review to positive_review

comment:79 Changed 3 years ago by rws

Thanks.

comment:80 Changed 3 years ago by vbraun

  • Branch changed from u/rws/16697-2 to 157b268d4ac45de049089e633809f89d2c3ab84e
  • Resolution set to fixed
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.