Opened 4 years ago

Closed 2 years ago

Last modified 18 months ago

#11143 closed defect (fixed)

define symbolic functions for exponential integrals

Reported by: kcrisman Owned by: benjaminfjones
Priority: major Milestone: sage-5.3
Component: symbolics Keywords: ei Ei special function maxima sd32 sd40.5
Cc: Merged in: sage-5.3.beta2
Authors: Benjamin Jones, Volker Braun Reviewers: Burcin Erocal, Karl-Dieter Crisman, William Stein
Report Upstream: N/A Work issues:
Branch: Commit:
Dependencies: #13109 Stopgaps:

Description (last modified by benjaminfjones)

We're missing some conversions from Maxima. Like exponential integrals of various kinds.

sage: f(x) = e^(-x) * log(x+1)
sage: uu = integral(f,x,0,oo)
sage: uu
x |--> e*expintegral_e(1, 1)

See this ask.sagemath post for some details.

Current symbol conversion table

From sage.symbolic.pynac.symbol_table['maxima'] as of Sage-4.7

Maxima ---> Sage

%gamma ---> euler_gamma
%pi ---> pi
(1+sqrt(5))/2 ---> golden_ratio
acos ---> arccos
acosh ---> arccosh
acot ---> arccot
acoth ---> arccoth
acsc ---> arccsc
acsch ---> arccsch
asec ---> arcsec
asech ---> arcsech
asin ---> arcsin
asinh ---> arcsinh
atan ---> arctan
atan2 ---> arctan2
atanh ---> arctanh
binomial ---> binomial
brun ---> brun
catalan ---> catalan
ceiling ---> ceil
cos ---> cos
delta ---> dirac_delta
elliptic_e ---> elliptic_e
elliptic_ec ---> elliptic_ec
elliptic_eu ---> elliptic_eu
elliptic_f ---> elliptic_f
elliptic_kc ---> elliptic_kc
elliptic_pi ---> elliptic_pi
exp ---> exp
expintegral_e ---> En
factorial ---> factorial
gamma_incomplete ---> gamma
glaisher ---> glaisher
imagpart ---> imag_part
inf ---> +Infinity
infinity ---> Infinity
khinchin ---> khinchin
kron_delta ---> kronecker_delta
li[2] ---> dilog
log ---> log
log(2) ---> log2
mertens ---> mertens
minf ---> -Infinity
psi[0] ---> psi
realpart ---> real_part
signum ---> sgn
sin ---> sin
twinprime ---> twinprime

Summary of missing conversions

Special functions defined in Maxima

(http://maxima.sourceforge.net/docs/manual/en/maxima_16.html#SEC56)

bessel_j (index, expr)         Bessel function, 1st kind
bessel_y (index, expr)         Bessel function, 2nd kind
bessel_i (index, expr)         Modified Bessel function, 1st kind
bessel_k (index, expr)         Modified Bessel function, 2nd kind
  • Notes: bessel_I, bessel_J, etc. are functions in Sage for numerical evaluation. There is also the Bessel class, but no conversions from Maxima's bessel_i etc. to Sage.
hankel_1 (v,z)                 Hankel function of the 1st kind
hankel_2 (v,z)                 Hankel function of the 2nd kind
struve_h (v,z)                 Struve H function
struve_l (v,z)                 Struve L function
  • Notes: None of these functions are currently exposed at the top level in Sage. Evaluation is possible using mpmath.
assoc_legendre_p[v,u] (z)      Legendre function of degree v and order u 
assoc_legendre_q[v,u] (z)      Legendre function, 2nd kind
  • Notes: In Sage we have legendre_P(n, x) and legendre_Q(n, x) both described as Legendre functions. It's not clear to me how there are related to Maxima's versions since the number of arguments differs.
%f[p,q] ([], [], expr)         Generalized Hypergeometric function
hypergeometric(l1, l2, z)      Hypergeometric function
slommel
%m[u,k] (z)                    Whittaker function, 1st kind
%w[u,k] (z)                    Whittaker function, 2nd kind
  • Notes: hypergeometric(l1, l2, z) needs a conversion to Sage's hypergeometric_U. The others can be evaluated using mpmath. slommel is presumably mpmath's lommels1() or lommels2() (or both?). This isn't well documented in Maxima.
expintegral_e (v,z)            Exponential integral E
expintegral_e1 (z)             Exponential integral E1
expintegral_ei (z)             Exponential integral Ei
expintegral_li (z)             Logarithmic integral Li
expintegral_si (z)             Exponential integral Si
expintegral_ci (z)             Exponential integral Ci
expintegral_shi (z)            Exponential integral Shi
expintegral_chi (z)            Exponential integral Chi
erfc (z)                       Complement of the erf function
  • Notes: The exponential integral functions expintegral_e1 and expintegral_ei (z) are called exponential_integral_1 and Ei resp. in Sage. They both need conversions. The rest need BuiltinFunction classes defined for them with evaluation handled by mpmath and the symbol table conversion added. Also, erfc is called error_fcn, so also needs a conversion.
kelliptic (z)                  Complete elliptic integral of the first 
                               kind (K)
parabolic_cylinder_d (v,z)     Parabolic cylinder D function
  • Notes: kelliptic(z) needs a conversion to elliptic_kc in Sage and parabolic_cylinder_d (v,z) does not seem to be exposed at top level. It can be evaluated by mpmath.

Apply to the Sage library:

  1. trac_11143-v2.5-rebased.4.patch
  2. trac-11143-ref.2.patch
  3. trac_11143-pynac-serials.patch

Attachments (15)

trac_11143_En.patch (7.3 KB) - added by benjaminfjones 3 years ago.
added the symbolic function exp_integral_e
trac_11143.patch (53.7 KB) - added by benjaminfjones 3 years ago.
adds sage/functions/exp_integral.py
trac_11143_v2.patch (50.7 KB) - added by benjaminfjones 3 years ago.
renamed functions according to trac comments and rebased to 5.0.beta5
trac_11143_v2.2.patch (50.7 KB) - added by benjaminfjones 3 years ago.
rebase to 5.0.beta10
trac_11143_v2.3.patch (50.7 KB) - added by benjaminfjones 3 years ago.
removed trailing whitespace
trac_11143_v2.4.patch (51.2 KB) - added by benjaminfjones 3 years ago.
trac_11143_v2.5.patch (51.5 KB) - added by benjaminfjones 3 years ago.
addresses log_integral(0) -> 0, exponential_integral_1(0) -> Infinity, sinh_integral issues, etc..
trac-11143-ref.2.patch (2.0 KB) - added by benjaminfjones 3 years ago.
trac-11143-ref.patch (2.0 KB) - added by was 3 years ago.
fixed
trac_11143-v2.5-rebased.patch (51.6 KB) - added by kcrisman 2 years ago.
Rebased to 5.1.beta6
trac_11143-pynac-serials.patch (4.5 KB) - added by kcrisman 2 years ago.
trac_11143-v2.5-rebased.2.2.patch (51.9 KB) - added by vbraun 2 years ago.
Rediffed patch
trac_11143-v2.5-rebased.2.patch (51.9 KB) - added by vbraun 2 years ago.
Rebased patch
trac_11143-v2.5-rebased.3.patch (51.8 KB) - added by benjaminfjones 2 years ago.
fixed failing doctest
trac_11143-v2.5-rebased.4.patch (51.8 KB) - added by benjaminfjones 2 years ago.
changed abs tol to rel tol

Download all attachments as: .zip

Change History (112)

comment:1 follow-ups: Changed 4 years ago by benjaminfjones

As far as I can tell, the general exponential_e function isn't available directly in Sage or in PARI (which is used to evaluate the exponential_integral_1 function in Sage).

Also, it's possible to get maxima to rewrite the exponential integrals in terms of gamma functions like so:

sage: maxima.eval('expintrep:gamma_incomplete')
'gamma_incomplete'
sage: maxima.integrate(exp(-x)*log(x+1), x, 0, oo)
%e*gamma_incomplete(0,1)
sage: N(e*gamma(0,1), digits=18)
0.596347362323194074

But as you see, gamma_incomplete isn't defined in Sage either, but the table sage.symbolic.pynac.symbol_table['maxima'] lists the Sage equivalent gamma. Anyway, it should be possible to have the maxima interface (with the help of maxima itself) rewrite any exponential integral that Sage doesn't have in terms of gamma functions.

By the way, the owner on the ticket is @burcin, does that mean they are working on it currently?

comment:2 Changed 4 years ago by benjaminfjones

  • Cc benjaminfjones added

comment:3 in reply to: ↑ 1 Changed 4 years ago by kcrisman

Replying to benjaminfjones:

As far as I can tell, the general exponential_e function isn't available directly in Sage or in PARI (which is used to evaluate the exponential_integral_1 function in Sage).

That's unfortunate. However, mpmath seems to have it. So we could create a symbolic version and have the _eval_ method call mpmath, which we seem to be moving to.

Also, it's possible to get maxima to rewrite the exponential integrals in terms of gamma functions like so:

sage: maxima.eval('expintrep:gamma_incomplete')
'gamma_incomplete'
sage: maxima.integrate(exp(-x)*log(x+1), x, 0, oo)
%e*gamma_incomplete(0,1)
sage: N(e*gamma(0,1), digits=18)
0.596347362323194074

Interesting.

But as you see, gamma_incomplete isn't defined in Sage either, but the table sage.symbolic.pynac.symbol_table['maxima'] lists the Sage equivalent gamma.

That should be ok; the whole point of the table is to convert into the Sage equivalent.

By the way, the owner on the ticket is @burcin, does that mean they are working on it currently?

No, that is an automatic thing that happens. It is possible to be designated an 'owner' of a ticket in a given component, which basically means you automatically get updates. If you want to 'own' it, please do! We really have plenty of special functions in Sage, but they are not always well exposed at the top level.

Incidentally, once you comment on a ticket, I believe the default is to copy you in on all replies. So you didn't have to cc: yourself specially :)

comment:4 in reply to: ↑ 1 Changed 4 years ago by burcin

Replying to benjaminfjones:

But as you see, gamma_incomplete isn't defined in Sage either, but the table sage.symbolic.pynac.symbol_table['maxima'] lists the Sage equivalent gamma. Anyway, it should be possible to have the maxima interface (with the help of maxima itself) rewrite any exponential integral that Sage doesn't have in terms of gamma functions.

Incomplete gamma is defined in Sage. You can access it directly though incomplete_gamma() or gamma_inc(). The top level function gamma() behaves like incomplete gamma if you give it two arguments. IIRC, this is similar to maple.

By the way, the owner on the ticket is @burcin, does that mean they are working on it currently?

I am not working on it. The ticket status assigned is supposed to indicate that the owner is working on the problem, but we don't use that much either.

comment:5 follow-up: Changed 4 years ago by burcin

I guess the point of this ticket is to define symbolic function in Sage to represent exponential integrals, etc. The symbolic function class handles adding stuff to the conversion table automatically.

Can we replace this ticket with several beginner tickets? One for each function we are missing.

comment:6 in reply to: ↑ 5 Changed 4 years ago by kcrisman

Replying to burcin:

I guess the point of this ticket is to define symbolic function in Sage to represent exponential integrals, etc. The symbolic function class handles adding stuff to the conversion table automatically.

Can we replace this ticket with several beginner tickets? One for each function we are missing.

Well, that would be nice. But we could also presumably do it directly, if that would solve this problem for now. Well, either is fine as long as it were to happen. If you do split this, be sure to give a good template (I mean a link to the template, not write it yourself).

comment:7 Changed 4 years ago by benjaminfjones

I'm attempting to write a template for the expintegral_e function (denoted E_n(z) in A&S). As I'm looking through the code, I see several models used for the functions and classes in sage/functions/special.py and sage/functions/transcendental.py

  • Functions like Function_exp_integral (also called Ei) are defined as classes that inherit from BuiltInFunction and call the mpmath implementation when evaluated. The function DickmanRho also does this and includes other nice methods for approximating values and power series.
  • Functions like EllipticE inherit from MaximaFunction which handles evaluation, etc. through Maxima. It seems there is an advantage to the mpmath implementation because presumably the interface is faster and the precision is arbitrary (whereas Maxima is limited to 53 bits).
  • Functions like Li and error_fcn are simply wrapper functions that try to evaluate the input symbolically or numerically depending on context.
  • In sage/functions/trig.py there is a mixture of classes that derive from GinacFunction (and include information in their __init__ methods about conversions to other systems like Maxima or Mathematica) and also functions that derive from BuiltinFunction. It's not clear to me why some functions are Ginac and some are Builtin.

Questions:

  • For the purposes of this ticket, what is recommended? @kcrisman 's comment leads me to believe that inheriting from the BuiltinFunction class and using mpmath for evaluation is preferable.
  • Where should the various exponential integral special functions that we are missing go? sage/functions/special.py, sage/functions/transcendental.py, or somewhere else?

comment:8 Changed 4 years ago by kcrisman

Those are good questions.

I think the most important thing is to make sure that whatever is implemented has

  • good numerical evaluation (perhaps via mpmath)
  • translates back and forth to Maxima properly (for integration and limits)

My sense is that BuiltinFunction would be best. GinacFunction probably is only good for things that in fact evaluate in Ginac. This explains trig.py. Ginac page shows that the ones which are GinacFunctions are exactly the ones which Ginac in fact has. I don't think it has most of these other functions.

As for MaximaFunction, it seems to inherit from BuiltinFunction and lives in the special functions file. This dates from the days when Sage had very few options for symbolic stuff and evaluation, and so it just does a few extra things. If we moved to mpmath for a given function, we would probably use BuiltinFunction and then add evaluation options for Maxima and add to the conversion dictionary as needed.

In fact, it wouldn't be a bad idea to have two different eval procedures if possible...


As for where such things go, probably it would make sense to separate a lot of these special functions out into a separate file. The distinction between transcendental and special is not totally obvious, for instance :)


By the way, we certainly want fewer of the wrapper functions! But that will take a lot of tedious work.

comment:9 Changed 4 years ago by kcrisman

By the way, adding this will really be a great help. Sage has so many of the things almost anybody needs, but if you have to use mpmath or GSL or R or something else directly, it sort of makes Sage a moot point. The idea is having a one-stop shop.

comment:10 Changed 4 years ago by benjaminfjones

I've uploaded a first shot at a template for the functions we want to add in this ticket. The patch exposes the general complex exponential integral function En also called Function_exp_integral_En by adding a class to sage/functions/special.py (I can change where it goes later on if needed / desired). Numerical evaluation is handled by mpmath and symbolics are handled by Sage (e.g. the derivative) and Maxima (e.g. the antiderivative).

One of the docstring examples shows that the integral of e^(-x) * log(x+1) from the ticket description is now evaluated properly.

Any comments or suggestions?

comment:11 Changed 4 years ago by benjaminfjones

  • Description modified (diff)
  • Keywords special function maxima added

comment:12 Changed 4 years ago by burcin

  • Authors set to Benjamin Jones
  • Milestone set to sage-4.7.1
  • Reviewers set to Burcin Erocal
  • Summary changed from Add various Maxima special functions to symbol table to define symbolic functions for exponential integrals

The patch looks great. Thanks for the template. A few suggestions:

  • The function name should be more explicit. I suggest exp_integral_e.
  • In the top level name space En can be an alias for this function (though I'd prefer not to take up a two letter name), but we should have the long name available. It is easier to find all these functions by exp_integral<tab> than E<tab>.
  • See _eval_ method of sage.functions.other.Function_gamma_inc for an example of how to find a common parent for the arguments
  • The call to mpmath does not require the prec parameter any more. If you pass a Sage element to mpmath it handles the precision correctly. You can delete the code for this in _evalf_().
  • You should remove the __call__() method altogether. The only purpose of that is to display the deprecation notice. Since you are implementing a new function here, there is nothing being deprecated.
  • In _derivative_() the call to this function should be exp_integral_e(n-1, z), assuming you change the function name accordingly.

I changed the ticket description to limit this to implementing symbolic functions for exponential integrals. We can use the wiki page for a general overview of the progress on symbolic functions.

comment:13 Changed 4 years ago by benjaminfjones

Thanks for the feedback, Burcin! I feel like I'm developing an understanding for how symbolic functions are handled in Sage now.

To your suggestions:

  • I agree, the changed the name of the class to Function_exp_integral_e and the global function name to exp_integral_e. I think it would make sense as part of this ticket to add an alias of exp_integral_1 for the existing exponential_integral_1 to be consistent.
  • Good point, I didn't think of that.
  • The function now correctly coerces the two arguments to the same parent. I also added automatic evaluation of some special cases (n=0 or z=0) as in the gamma_inc function.
  • I played around with removing the prec argument, but didn't get the results I expected. I found that passing the parent explicitly works. The common parent is determined in _eval_ (as in the gamma_inc function).
  • Done.
  • Done.

Let me know what you think. Once we've agreed on a good template for this one function, I will implement the rest of the exponential integral functions that are now listed on the wiki and upload to this ticket.

comment:14 Changed 4 years ago by burcin

Perfect! I can only nitpick about documentation. :)

  • lines should be < 80 characters long. This helps when viewing help in a terminal.
  • I don't know how testing for 0 in _eval_ effects performance. This is a hard problem. See the _eval_ methods in sage/functions/generalized.py for a workaround if this is too slow.

You're right, we should add an alias exp_integral_1.

Thanks for working on this.

comment:15 follow-up: Changed 4 years ago by benjaminfjones

  • Cc benjaminfjones removed

The tests similar to if z == 0 in _eval_ do make a big difference. I guess this is well known, but I didn't realize how big the speed difference is.

I made a little table below of some timings. In the first table , _eval_ includes tests if z == 0 and n > 1 and if n == 0. In the second table, there are no such if statements (so the special cases are not implemented). In the third table, the special cases are implementes as in sage/functions/generalized.py where an approximation of z (or n) is produced and checked instead of the symbol.

with ``if z == 0``

=============================================   ========
  Test                                                                                   Time
=============================================   ========
sage: timeit("f = exp_integral_e(n,z)")                  1.44 ms
sage: timeit("f = exp_integral_e(n,0)")                  929 µs
sage: timeit("f = exp_integral_e(0,z)")                  1.41 ms
sage: timeit("f = exp_integral_e(1.0,1.0)")              158 µs
=============================================   ========

without ``if z == 0``

=============================================   ======
  Test                                                                                   Time
=============================================   ======
sage: timeit("f = exp_integral_e(n,z)")                 541 µs
sage: timeit("f = exp_integral_e(n,0)")                 300 µs
sage: timeit("f = exp_integral_e(0,z)")                 299 µs
sage: timeit("f = exp_integral_e(1.0,1.0)")             161 µs
=============================================   ======

with:

.. code-block:: python

        try:
            approx_z = ComplexIntervalField()(z)
            # if z is zero and n > 1
            if bool(approx_z.imag() == 0) and bool(approx_z.real() == 0):
                if n > 1:
                    return 1/(n-1)
        except: # z is symbolic
            pass
        # if n is zero
        try:
            approx_n = ComplexIntervalField()(n)
            if bool(approx_n.imag() == 0) and bool(approx_n.real() == 0):
                return exp(-z)/z
        except: # n is symbolic
            pass

=============================================   ======
  Test                                                                                   Time
=============================================   ======
sage: timeit("f = exp_integral_e(n,z)")                 570 µs
sage: timeit("f = exp_integral_e(n,0)")                 576 µs
sage: timeit("f = exp_integral_e(0,z)")                 1.05 ms
sage: timeit("f = exp_integral_e(1.0,1.0)")             160 µs
=============================================   ======

Timings in tables 2 and 3 are close except in the case where exp(-z)/z is
returned, whereas table 1 is anywhere from a factor of 3 to 5 slower than in table 2 when a symbolic argument is passed. Anyway, I thought I'd include the above for other beginners such as myself.


Another thing I discovered is that these two special cases that I was
implementing are known to Maxima:

sage: f = exp_integral_e(0,x)
sage: f
exp_integral_e(0,x)
sage: f.simplify()
e^(-x)/x

sage: nn = var('nn')
sage: assume(nn > 1)
sage: f = exp_integral_e(nn, 0)
sage: f
exp_integral_e(nn, 0)
sage: f.simplify()
1/(nn - 1)


So I think in the interest of speedy evaluation it's best to leave the special
cases out, but point out in the documentation that Maxima knows about them.

I've uploaded a new patch. I'll start implementing the other exponential integrals using this as a template.

comment:16 in reply to: ↑ 15 Changed 4 years ago by burcin

The timings are really instructive, we should link to them from wiki page. Thank you for this careful study.

Replying to benjaminfjones:

So I think in the interest of speedy evaluation it's best to leave the special
cases out, but point out in the documentation that Maxima knows about them.

I think your timings show that using the CIF approximation is not such a big penalty. Note that your timings also reflect the construction of the result and not just including that change. It is better if Sage can do the evaluation automatically without having to pass through a simplify() call. I'd like them in, with the approximation approach.

Since this approximation is being called so often, we should make it a method of symbolic expressions. I opened a new ticket for this: #11513. This ticket should depend on that.

I will provide a preliminary patch for #11513 soon.

comment:17 Changed 4 years ago by burcin

I attached a very simple patch to #11513. Now we can write the check for zero as:

def is_zero(z):
    if isinstance(z, Expression):
        return z._is_numerically_zero()
    else:
        return not z

It would be nice if we could simplify this further, but I need to go and do other things now. :)

comment:18 Changed 3 years ago by benjaminfjones

The patch now depends on #11513

Here are new timings for the _eval_ method with _is_numerically_zero():

.. code-block:: python

    # special case: *quickly* test if (z == 0 and n > 1)
    if isinstance(z, Expression):
        if z._is_numerically_zero():
            z_zero = True # for later
            if n > 1:
                return 1/(n-1)
        else:
            if not z: 
                z_zero = True
                if n > 1:
                    return 1/(n-1)

======================================================  =======
  Test                                                   Time
======================================================  =======
sage: timeit("f = exp_integral_e(n,z)")                 535 µs
sage: timeit("f = exp_integral_e(n,0)")                 482 µs
sage: assume(n > 1); timeit("f = exp_integral_e(n,0)")  3.56 ms
sage: timeit("f = exp_integral_e(0,z)")                 968 µs
sage: timeit("f = exp_integral_e(1.0,1.0)")             160 µs
======================================================  =======

I realized that in row 2 of the previous timings I neglected to assume n > 1 so those timings aren't giving much information since the expression is left unevaluated like in row 1. The new row 3 includes that assumption so that the simplified result 1/(n-1) is created and returned.

I'll update the timings above and move these tables to the wiki.

Changed 3 years ago by benjaminfjones

added the symbolic function exp_integral_e

comment:19 follow-up: Changed 3 years ago by benjaminfjones

  • Owner changed from burcin to benjaminfjones

In the patch, I'm added a function to sage/functions/special.py. But with the six new functions, this is going to be *a lot* of code in the one file.

I was thinking of moving the six new functions to a new file sage/functions/exp_integral.py. Any thoughts? If that seems like a good idea, does it also make sense to move the exp integrals that already exist in Sage to the same place (e.g. exponential_integral_1 lives now in sage/functions/transcendental.py.

comment:20 in reply to: ↑ 19 Changed 3 years ago by kcrisman

I was thinking of moving the six new functions to a new file sage/functions/exp_integral.py. Any thoughts? If that seems like a good idea, does it also make sense to move the exp integrals that already exist in Sage to the same place (e.g. exponential_integral_1 lives now in sage/functions/transcendental.py.

Yes, absolutely. The functions/ folder needs to be reorganized, though obviously this is lower priority. Make sure that it still shows up in the reference manual, that all imports work fine, etc. But these functions are all supposed to be accessed through the top level namespace anyway, so this is a great idea to make it easier to organize.

comment:21 Changed 3 years ago by benjaminfjones

I've uploaded a new partial patch for the ticket. I wanted to post this now before I finish the patch so people can comment on the consolidation of exponential functions in the new file sage/functions/exp_integral.py and/or anything else.

If the basic organisation looks good, I'll finish up adding the remaining exponential integrals and upload a new patch.

comment:22 Changed 3 years ago by kcrisman

  • Dependencies set to #11513
  • Reviewers changed from Burcin Erocal to Burcin Erocal, Karl-Dieter Crisman

A few comment:

  • This will conflict (slightly) with #11043, I think. You can help their cause by doing
    from sage.plot.all import plot
    
    I believe. The point is making sure we don't import numpy and other things like that on startup. It's possible that some of the other imports will have this problem, but hopefully not (?).
  • I don't really understand coercion, so I'm reluctant to comment on the coercion stuff in _eval_, though I'd be glad if you explained that (I don't recall seeing it in these types of functions before).
  • This seems to depend on #11513 still, so I'm adding that. But that still does not have positive review. Just pointing it out.
  • Line 260: raise NotImplementedError("The derivative of is only implemented for n = 1, 2, 3, ...") seems to be missing something.

On a cursory reading, this looks really nice. Certainly you should go ahead with the full reorganization, that will help a lot.

comment:23 Changed 3 years ago by benjaminfjones

  • Keywords sd32 added

comment:24 follow-up: Changed 3 years ago by benjaminfjones

Trying to finish up the patch, I've run into a problem evaluating any of my new functions at python floats. For example, apply the patch and try:

sage: exp_integral_e1(float(1))
ERROR: An unexpected error occurred while tokenizing input
The following traceback may be corrupted or invalid
The error message is: ('EOF in multi-line statement', (690, 0))

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)

/Users/jonesbe/sage/sage-4.7.2.alpha2/devel/sage-test/sage/<ipython console> in <module>()

/Users/jonesbe/sage/latest/local/lib/python2.6/site-packages/sage/symbolic/function.so in sage.symbolic.function.Function.__call__ (sage/symbolic/function.cpp:4599)()

/Users/jonesbe/sage/latest/local/lib/python2.6/site-packages/sage/functions/exp_integral.pyc in _eval_(self, z)
    335         """
    336         if not isinstance(z, Expression) and is_inexact(z):
--> 337             return self._evalf_(z, parent(z))
    338 
    339         return None # leaves the expression unevaluated

/Users/jonesbe/sage/latest/local/lib/python2.6/site-packages/sage/functions/exp_integral.pyc in _evalf_(self, z, parent)
    350         """
    351         import mpmath
--> 352         return mpmath_utils.call(mpmath.e1, z, parent=parent)
    353 
    354     def _derivative_(self, z, diff_param=None):

/Users/jonesbe/sage/latest/local/lib/python2.6/site-packages/sage/libs/mpmath/utils.so in sage.libs.mpmath.utils.call (sage/libs/mpmath/utils.c:5277)()

AttributeError: type object 'float' has no attribute 'prec'

You also get the error if you try to plot any of the exponential integral functions as they are written unless you wrap the input with RDF or something similar.

The error is being raised inside the call function in sage/libs/mpmath/util.pyx. I've looked through that code and thought about what would happen if the prec option is not specified (as in my code) and the parent is None (as when a python float is passed), but I can't figure out why it doesn't work to pass a float.

On the other hand, changing the lines in the _evalf method

import mpmath
return mpmath_utils.call(mpmath.e1, z, parent=parent)

to

import mpmath
if isinstance(parent, Parent) and hasattr(parent, 'prec'):
    prec = parent.prec()
else:
    prec = 53
return mpmath_utils.call(mpmath.si, z, prec=prec)

fixes the error.

What do people think? Is this a bug in the call function (which claims it will accept any type that mpmath handles natively like python floats) or is there something wrong with my code?

comment:25 in reply to: ↑ 24 Changed 3 years ago by burcin

Replying to benjaminfjones:

What do people think? Is this a bug in the call function (which claims it will accept any type that mpmath handles natively like python floats) or is there something wrong with my code?

It would be better to fix the mpmath_call() function. This is called from many places and duplicating the code to check if parent has a prec() method everywhere is error prone.

Thanks for your work on this so far. I will try to come up with an acceptable implementation for #11513 next week. I hope it doesn't hold this up for long.

comment:26 Changed 3 years ago by benjaminfjones

Ok, I've created a new ticket #11885 to fix the issue with mpmath.utils.call and uploaded a patch there. I'll upload the finished patch for #11143 that depends on it later this afternoon.

comment:27 Changed 3 years ago by benjaminfjones

  • Description modified (diff)

comment:28 Changed 3 years ago by benjaminfjones

  • Status changed from new to needs_review

The patch just uploaded should be feature complete. I've checked the documentation, checked that the doctests pass, and checked that make ptestlong completes without errors on Mac OS X 10.6.8 and also debian linux squeeze. The patch is based on sage-4.7.2.alpha3.

Changed 3 years ago by benjaminfjones

adds sage/functions/exp_integral.py

comment:29 Changed 3 years ago by benjaminfjones

There were a few doctests failing in sage/misc/sagedoc.py that I missed. The patch just uploaded now *really* doesn't break anything according to make ptestlong.

There are some strange doctests in sagedoc.py, especially the recursive ones. I fixed doctests that need changing when a pynac function is added to the Sage library and then, because I changed a docstring in that file, another doctest failed. Phew..

comment:30 Changed 3 years ago by kcrisman

  • Dependencies changed from #11513 to #11513, 11885

Impressive work. No guarantees that I'll be able to look through this long patch with a fine-tooth comb anytime soon, but this will be a great addition - I really like Ci and Si being in.

    \operatorname{li}(x) = \int_0^z \frac{dt}{\ln(t)} = \operatorname{Ei}(\ln(x)) 

?? What happened to the x in the middle?

Also, what relation does this function bear to Li (I assume it's the \int_0 as opposed to \int_2, if I recall correctly), and why is that not mentioned? See #7357. I have to say that the names of the functions are ... very long. Is there a problem with the fairly standard Li and li?

See also #3401, though I don't know if that is directly relevant for this particular set of functions.

comment:31 Changed 3 years ago by kcrisman

  • Dependencies changed from #11513, 11885 to #11513, #11885

comment:32 Changed 3 years ago by benjaminfjones

@kcrisman - Yes, thanks, I just spotted that typo in the docstring for exp_integral_li while I was commenting on the other tickets you cc'd me on.

The function names are long. I wanted to choose a common prefix for all the exponential integrals. I think burcin suggested exp_integral originally. I think the problem with li and Li, or En and Ei, etc.. is that the name is very short and not very informative unless you are familiar with all the different varieties of standard exponential integrals. Also, I think the short two-letter names would be easy to mix up, especially li and Li. Yes, you're right, the difference between those logarithmic integrals is the constant \int_0^2 1/log(t) dt.

I made some comments on #7357 and #3401. Basically, I think it makes sense to roll the issues in those two tickets into this one since this one is an overhaul of the exponential integral functions with the aim to make them all symbolic and organize them in the same module.

comment:33 Changed 3 years ago by benjaminfjones

Sorry for the double comment. I remembered the rational for the common prefix. Burcin point out that once a user knows about one exponential integral function, they can type:exp_integral<tab> and see the list of all of them. That wouldn't be possible if the names are just the standard mathematical notation li, Li, En, Ei, etc..

comment:34 Changed 3 years ago by kcrisman

That second comment about tab-completion makes sense. I just think it would be helpful to have easy-to-find aliases as well. Here are Mathematica's names. I understand the concern about the standard notation, but it is standard, and I would never think to look for log integral under exp_integral<tab>, if you catch my drift - and I do know it can be written as one.

Sorry about the comments on the other tickets about naming conventions - feel free to ignore those, I read this update last.

comment:35 follow-up: Changed 3 years ago by benjaminfjones

I see your point. How about the following names:

  • exp_integral_e
  • exp_integral_e1
  • exp_integral_ei
  • log_integral (for li)
  • log_integral_offset for (Li)
  • cos_integral
  • sin_integral
  • sinh_integral
  • cosh_integral

... and maybe aliases for the usual mathematical names: li, Li, En, Ei, E1, si, ci, shi, chi ? Perhaps that would be too many new names in the global name space?

comment:36 in reply to: ↑ 35 Changed 3 years ago by kcrisman

  • Status changed from needs_review to needs_work

Replying to benjaminfjones:

I see your point. How about the following names:

  • exp_integral_e
  • exp_integral_e1
  • exp_integral_ei
  • log_integral (for li)
  • log_integral_offset for (Li)
  • cos_integral
  • sin_integral
  • sinh_integral
  • cosh_integral

This could work, and is easier to understand. Maybe exp_integral<tab> could also bring up some command which gave general help for this whole module... ?

... and maybe aliases for the usual mathematical names: li, Li, En, Ei, E1, si, ci, shi, chi ? Perhaps that would be too many new names in the global name space?

You're right that it's a lot. On the other hand, Si and li/Li are fairly standard... I think it might be worth asking on sage-devel as to how many of these would be useful. I definitely think that things which are global in both Mma and Maple should be okay. See Maple help and e.g. this Mathematica help. Seems like Maple is more likely to support the standard notation as a shortcut. But at the very least your ideas above are better than the initial one.

I'll put this as 'needs work' for this reason.


Thank you so much for the hard work and at times surely boring work improving this. There are so many things already in Sage which we just need to make easier for the end user.

Changed 3 years ago by benjaminfjones

renamed functions according to trac comments and rebased to 5.0.beta5

comment:37 follow-up: Changed 3 years ago by benjaminfjones

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

Since all the dependencies have been merged I thought I would rebase this ticket and make the outstanding changes from the last comment. I've changed the names of the functions in the new module sage/functions/exp_integral.py according to my proposal in comment 36 and added top-level aliases li, si, ci, shi, chi.

comment:38 in reply to: ↑ 37 Changed 3 years ago by kcrisman

Since all the dependencies have been merged I thought I would rebase this ticket and make the outstanding changes from the last comment. I've changed the names of the functions in the new module sage/functions/exp_integral.py according to my proposal in comment 36 and added top-level aliases li, si, ci, shi, chi.

This came up while preparing for a Sage talk today - I tried

sage: integrate(log(x-1)*e^x,x)
e*expintegral_e(1, -x + 1) + e^x*log(x - 1)

and of course it doesn't LaTeX very nicely.

I'll try to look at this soon. I already know that we would need to make si into Si, etc., because that is the standard notation (first letter capitalized, see all the references above), but I can make a reviewer patch which does this. But it sounds like this should be finally ready to go. Sorry for the long process here!

comment:39 Changed 3 years ago by benjaminfjones

I can fix the latex issues if you haven't started already. I have the afternoon off.

comment:40 Changed 3 years ago by kcrisman

I don't know that there are any LaTeX issues; the example above was with Sage 4.8! I suppose you can try that example with this patch to see.

I only commented to remind myself that we should have Si and not si, in line with Maple and with Mathematica's documentation on the 'usual' name for this function (though not its command). It keeps this patch on my mind. If there is something wrong, though, you can definitely fix it!

comment:41 Changed 3 years ago by davidloeffler

Apply trac_11143_v2.patch

(for patchbot)

comment:42 follow-up: Changed 3 years ago by davidloeffler

  • Status changed from needs_review to needs_work
  • Work issues set to rebase to 5.0.beta10

This patch doesn't apply to the current Sage beta -- see patchbot logs.

Changed 3 years ago by benjaminfjones

rebase to 5.0.beta10

comment:43 in reply to: ↑ 42 Changed 3 years ago by benjaminfjones

Replying to davidloeffler:

This patch doesn't apply to the current Sage beta -- see patchbot logs.

Thanks for updating the ticket description. I've rebased my patch to 5.0.beta10.

comment:44 Changed 3 years ago by davidloeffler

  • Status changed from needs_work to needs_review

comment:45 Changed 3 years ago by JoalHeagney

Hi guys. I don't know if this even applies, because I discovered this behaviour on version 4.8, with the patch applied (And it seems to apply to limit rather than anything else.)
However, if I enter the following code:

def fracintegral(func,n, xsub=x,a=0):
    t = var('t')
    assume(t > a)
    assume(x > a)
    return integrate(((x-t)^(n-1))*func(t),t,a,x)/gamma(n)
a = fracintegral(sin(x),1/2);a

I get the following returned:

-1/2*(sqrt(x)*sin(x)*exp_integral_e(1/2, -I*x) +
sqrt(x)*sin(x)*exp_integral_e(1/2, I*x) +
I*sqrt(x)*cos(x)*exp_integral_e(1/2, -I*x) -
I*sqrt(x)*cos(x)*exp_integral_e(1/2, I*x) -
2*limit(1/2*sqrt(-t)*sin(x)*exp_integral_e(1/2, -I*t) +
1/2*sqrt(-t)*sin(x)*exp_integral_e(1/2, I*t) -
1/2*I*sqrt(-t)*cos(x)*exp_integral_e(1/2, -I*t) +
1/2*I*sqrt(-t)*cos(x)*exp_integral_e(1/2, I*t), t, 0, minus))/sqrt(pi)

If I copy/paste this into a new sage notebook cell, I get the following error:

Traceback (most recent call last):    1/2*I*sqrt(-t)*cos(x)*exp_integral_e(1/2, -I*t) +
  File "", line 1, in <module>
    
  File "/tmp/tmpbIxhPU/___code___.py", line 10, in <module>
    _sage_const_1 /_sage_const_2 *I*sqrt(-t)*cos(x)*exp_integral_e(_sage_const_1 /_sage_const_2 , I*t), t, _sage_const_0 , minus))/sqrt(pi)
NameError: name 'minus' is not defined

If I fix up the minus by doing this,

-1/2*(sqrt(x)*sin(x)*exp_integral_e(1/2, -I*x) +
sqrt(x)*sin(x)*exp_integral_e(1/2, I*x) +
I*sqrt(x)*cos(x)*exp_integral_e(1/2, -I*x) -
I*sqrt(x)*cos(x)*exp_integral_e(1/2, I*x) -
2*limit(1/2*sqrt(-t)*sin(x)*exp_integral_e(1/2, -I*t) +
1/2*sqrt(-t)*sin(x)*exp_integral_e(1/2, I*t) -
1/2*I*sqrt(-t)*cos(x)*exp_integral_e(1/2, -I*t) +
1/2*I*sqrt(-t)*cos(x)*exp_integral_e(1/2, I*t), t, 0, 'minus'))/sqrt(pi)

I get this error:

Traceback (most recent call last):    1/2*I*sqrt(-t)*cos(x)*exp_integral_e(1/2, -I*t) +
  File "", line 1, in <module>
    
  File "/tmp/tmpI2jAGB/___code___.py", line 10, in <module>
    _sage_const_1 /_sage_const_2 *I*sqrt(-t)*cos(x)*exp_integral_e(_sage_const_1 /_sage_const_2 , I*t), t, _sage_const_0 , 'minus'))/sqrt(pi)
  File "/home/joal/bin/sage-4.8-linux-64bit-ubuntu_10.04.3_lts-x86_64-Linux/local/lib/python2.6/site-packages/sage/calculus/calculus.py", line 1131, in limit
    raise ValueError, "call the limit function like this, e.g. limit(expr, x=2)."
ValueError: call the limit function like this, e.g. limit(expr, x=2).

Then if I write it up in correct limit syntax like so:

b = 2*limit(1/2*sqrt(-t)*sin(x)*exp_integral_e(1/2, -I*t) +
1/2*sqrt(-t)*sin(x)*exp_integral_e(1/2, I*t) -
1/2*I*sqrt(-t)*cos(x)*exp_integral_e(1/2, -I*t) +
1/2*I*sqrt(-t)*cos(x)*exp_integral_e(1/2, I*t), t=0, dir='minus'); b

I get literally nothing from the notebook cell. If I then enter b and evaluate in the next cell I get this:

2*limit(1/2*I*sqrt(t)*sin(x)*exp_integral_e(1/2, -I*t) +
1/2*I*sqrt(t)*sin(x)*exp_integral_e(1/2, I*t) +
1/2*sqrt(t)*cos(x)*exp_integral_e(1/2, -I*t) -
1/2*sqrt(t)*cos(x)*exp_integral_e(1/2, I*t), t, 0, minus)

This is all with Pretty Printing turned off. If I turn on Pretty Printing, then the reply from the fracintegral call (on trig functions and exponents) is this:

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "_sage_input_35.py", line 10, in <module>
    exec compile(u'open("___code___.py","w").write("# -*- coding: utf-8 -*-\\n" + _support_.preparse_worksheet_cell(base64.b64decode("YSA9IGZyYWNpbnRlZ3JhbChzaW4oeCksMS8yKTth"),globals())+"\\n"); execfile(os.path.abspath("___code___.py"))
  File "", line 1, in <module>
    
  File "/tmp/tmpOvxXQw/___code___.py", line 3, in <module>
    exec compile(u'a = fracintegral(sin(x),_sage_const_1 /_sage_const_2 );a
  File "", line 1, in <module>
    
  File "/home/joal/bin/sage-4.8-linux-64bit-ubuntu_10.04.3_lts-x86_64-Linux/local/lib/python2.6/site-packages/sage/misc/latex.py", line 2243, in pretty_print
    view(object)
  File "/home/joal/bin/sage-4.8-linux-64bit-ubuntu_10.04.3_lts-x86_64-Linux/local/lib/python2.6/site-packages/sage/misc/latex.py", line 1969, in view
    s = _latex_file_(objects, title=title, sep=sep, tiny=tiny, debug=debug, **latex_options)
  File "/home/joal/bin/sage-4.8-linux-64bit-ubuntu_10.04.3_lts-x86_64-Linux/local/lib/python2.6/site-packages/sage/misc/latex.py", line 1624, in _latex_file_
    L = latex(x)
  File "/home/joal/bin/sage-4.8-linux-64bit-ubuntu_10.04.3_lts-x86_64-Linux/local/lib/python2.6/site-packages/sage/misc/latex.py", line 875, in __call__
    return LatexExpr(x._latex_())
  File "expression.pyx", line 667, in sage.symbolic.expression.Expression._latex_ (sage/symbolic/expression.cpp:4121)
  File "ring.pyx", line 596, in sage.symbolic.ring.SymbolicRing._latex_element_ (sage/symbolic/ring.cpp:6461)
  File "pynac.pyx", line 433, in sage.symbolic.pynac.py_latex_function (sage/symbolic/pynac.cpp:4702)
  File "pynac.pyx", line 407, in sage.symbolic.pynac.py_latex_function_pystring (sage/symbolic/pynac.cpp:4328)
TypeError: _limit_latex_() takes exactly 4 arguments (5 given)

Let me know if I should put this into another ticket for limit, but kcrisman asked me to put this up here first.

Changed 3 years ago by benjaminfjones

removed trailing whitespace

comment:46 Changed 3 years ago by benjaminfjones

  • Dependencies #11513, #11885 deleted
  • Description modified (diff)

comment:47 Changed 3 years ago by zimmerma

  • Status changed from needs_review to needs_work
  • Work issues changed from rebase to 5.0.beta10 to define expintegral_ei

while looking in the sources for the old name expintegral_e I found the following
instance:

sage: search_src("expintegral_e")
functions/exp_integral.py:393:        x*log_integral(x) - expintegral_ei(2*log(x))
sage: expintegral_ei?
Object `expintegral_ei` not found.

That should be fixed too.

Paul

PS: the patch applied cleanly to 5.0, thus I removed the "rebase to 5.0.beta10" work issue

Changed 3 years ago by benjaminfjones

comment:48 Changed 3 years ago by benjaminfjones

  • Status changed from needs_work to needs_review
  • Work issues define expintegral_ei deleted

Good catch, I've added the maxima conversion for expintegral_ei in the new patch. I've also changed the top level aliases si -> Si, ci -> Ci, etc.. to match standard notation. The added file exp_integral.py passes all tests, I'm running the full test suite now.

Patchbot: apply trac_11143_v2.4.patch

Last edited 3 years ago by benjaminfjones (previous) (diff)

comment:49 Changed 3 years ago by benjaminfjones

  • Keywords sd40.5 added

comment:50 Changed 3 years ago by dsm

Looks good! A few quirks I've noticed:

(1) li behaves a little differently with respect to autoevaluation at 0:

sage: Si(0)
0
sage: Shi(0)
0
sage: li(0)
log_integral(0)
sage: li(0).simplify()
0
sage: li(0).n()
0.000000000000000

(2) In the (probably semi-deprecated; it works with Python floats) exponential_integral_1, we should special-case 0:

sage: exp_integral_e1(0.01)
4.03792957653811
sage: exponential_integral_1(0.01)
4.037929576538114
sage: exp_integral_e1(0)
exp_integral_e1(0)
sage: exp_integral_e1(0).n()
+infinity
sage: exponential_integral_1(0)
---------------------------------------------------------------------------
PariError                                 Traceback (most recent call last)

/Users/mcneil/sagedev/sage-5.1.beta0/devel/sage-hack/sage/<ipython console> in <module>()

/Users/mcneil/sagedev/sage-5.1.beta0/local/lib/python2.7/site-packages/sage/functions/exp_integral.pyc in exponential_integral_1(x, n)
   1327     from sage.libs.pari.all import pari
   1328     if n <= 0:
-> 1329         return float(pari(x).eint1())
   1330     else:
   1331         return [float(z) for z in pari(x).eint1(n)]

/Users/mcneil/sagedev/sage-5.1.beta0/local/lib/python2.7/site-packages/sage/libs/pari/gen.so in sage.libs.pari.gen._pari_trap (sage/libs/pari/gen.c:49907)()

PariError:  (5)

(3) There's something going on with integration I don't understand:

sage: integrate(sinh_integral(x), x)
x*sinh_integral(x) - cosh(x)
sage: integrate(sinh_integral(x), x, 0, 1/2)
-cosh(1/2) + 1/2*sinh_integral(1/2) + 1
sage: integrate(sinh_integral(x), x, 0, 1/2).n()
0.125872409703453
sage: 
sage: integrate(sinh_integral(x), x, 0.0, 0.5)
0.125872409703 + 1.57079632679*I

IOW, somehow we pick up a pi/2 I, and it probably happens on our side:

In [7]: quad(lambda x: shi(x), [0, 0.5])
Out[7]: mpf('0.12587240970345281')

comment:51 Changed 3 years ago by kcrisman

This is really great. I've checked a number of things against Wolfram Alpha, doc has lots of great examples, excellent work. Great catches, Doug.

We can now completely remove exp_int, because it was deprecated over two years ago.

See #3401 for making the offset log integral symbolic as well. William would really like this.

comment:52 follow-up: Changed 3 years ago by kcrisman

Okay, I've reported the issue with Maxima using the wrong branch for Shi at this bug tracker ticket.

(%i30) float(expintegral_shi(1/2));
(%o30) 3.141592653589793*%i+0.506996749819667

comment:53 Changed 3 years ago by dsm

Thanks, kcrisman. Since after real-world discussions it appears that the integration difficulty is another Maxima quirk, I don't think we should hold up approving this, especially because it can be worked around. I do think it's worth warning the user, and the example will serve as a reminder to us that it can be removed when we update to a version of Maxima which finally fixes it, 'cause the doctest will fail.

So I propose adding the following to EXAMPLES for sinh_integral:

    Note that due to some problems with the way Maxima handles these
    expressions, definite integrals can sometimes give unexpected
    results (typically when using inexact endpoints) due to
    inconsistent branching::

        sage: integrate(sinh_integral(x), x, 0, 1/2)
        -cosh(1/2) + 1/2*sinh_integral(1/2) + 1
        sage: integrate(sinh_integral(x), x, 0, 1/2).n() # right
        0.125872409703453
        sage: integrate(sinh_integral(x), x, 0, 0.5).n() # wrong!
        0.125872409703453 + 1.57079632679490*I

Changed 3 years ago by benjaminfjones

addresses log_integral(0) -> 0, exponential_integral_1(0) -> Infinity, sinh_integral issues, etc..

comment:54 Changed 3 years ago by benjaminfjones

The new patch should take care of dsm's (1), (2), (3) and kcrisman's exp_int removal. Please test the new automatic simplification code in log_integral and exponential_integral_1.

Patchbot: apply trac_11143_v2.5.patch

Last edited 3 years ago by benjaminfjones (previous) (diff)

comment:55 Changed 3 years ago by kcrisman

I don't like this.

sage: li(0)
0
sage: type(li(0))
<type 'sage.rings.integer.Integer'>
sage: Si(0)
0
sage: type(Si(0))
<type 'sage.rings.integer.Integer'>

The reason is

        # special case: z = 0 
	        if isinstance(z, Expression): 
	            if z.is_trivial_zero(): 
	                return z 
	        else: 
	            if not z: 
	                return z 

type code, because of course the zero z we input is not an expression, but is not z, so we return the Integer zero instead of the symbolic zero.

Of course, we're inconsistent here:

sage: sin(0)
0
sage: type(sin(0))
<type 'int'>

Oops. But under the philosophy of

sage: f(x) = x^2
sage: f(0)
0
sage: type(f(0))
<type 'sage.symbolic.expression.Expression'>

I think we should return an Expression.

comment:56 Changed 3 years ago by benjaminfjones

For the record, I think we should try as much as possible to return outputs with the same parent as the input (whenever this makes sense). For li(int(0)) we should get int(0), for li(Integer(0)) we should get Integer(0), etc..

comment:57 Changed 3 years ago by kcrisman

  • Description modified (diff)

comment:58 Changed 3 years ago by was

  • Reviewers changed from Burcin Erocal, Karl-Dieter Crisman to Burcin Erocal, Karl-Dieter Crisman, William Stein
  • Status changed from needs_review to positive_review

Looks good to me. I've attached a little referee patch.

comment:59 Changed 3 years ago by benjaminfjones

  • Description modified (diff)

Thanks! I'll update the ticket description.

comment:60 Changed 3 years ago by kcrisman

The reviewer patch introduces an erroneous hunk with a backtick in a random line.

Changed 3 years ago by benjaminfjones

comment:61 Changed 3 years ago by benjaminfjones

  • Description modified (diff)

Changed 3 years ago by was

fixed

comment:62 follow-up: Changed 3 years ago by JoalHeagney

With the patches applied to 5.0, my fractional integral code works on sine functions, but I still get the following error if pretty printing is on. Is there a trac post for this?

def fracintegral(func,xsub,n,a=0):
    t = var('t')
    assume(t > a)
    assume(x > a)
    return integrate((x-t)^(n-1)*func.subs({xsub:t}),t,a,x)/gamma(n)
a = fracintegral(sin(x),x,1/2)
a

Returns this on the call to a (but not on the assignment of a =):

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "_sage_input_19.py", line 10, in <module>
    exec compile(u'open("___code___.py","w").write("# -*- coding: utf-8 -*-\\n" + _support_.preparse_worksheet_cell(base64.b64decode("YQ=="),globals())+"\\n"); execfile(os.path.abspath("___code___.py"))
  File "", line 1, in <module>
    
  File "/tmp/tmpJfuEaN/___code___.py", line 2, in <module>
    exec compile(u'a
  File "", line 1, in <module>
    
  File "/home/joal/bin/sage-5.0/local/lib/python2.7/site-packages/sage/misc/latex.py", line 2280, in pretty_print
    view(object)
  File "/home/joal/bin/sage-5.0/local/lib/python2.7/site-packages/sage/misc/latex.py", line 2006, in view
    s = _latex_file_(objects, title=title, sep=sep, tiny=tiny, debug=debug, **latex_options)
  File "/home/joal/bin/sage-5.0/local/lib/python2.7/site-packages/sage/misc/latex.py", line 1661, in _latex_file_
    L = latex(x)
  File "/home/joal/bin/sage-5.0/local/lib/python2.7/site-packages/sage/misc/latex.py", line 909, in __call__
    return LatexExpr(x._latex_())
  File "expression.pyx", line 815, in sage.symbolic.expression.Expression._latex_ (sage/symbolic/expression.cpp:4898)
  File "ring.pyx", line 605, in sage.symbolic.ring.SymbolicRing._latex_element_ (sage/symbolic/ring.cpp:6658)
  File "pynac.pyx", line 433, in sage.symbolic.pynac.py_latex_function (sage/symbolic/pynac.cpp:4817)
  File "pynac.pyx", line 407, in sage.symbolic.pynac.py_latex_function_pystring (sage/symbolic/pynac.cpp:4443)
TypeError: _limit_latex_() takes exactly 4 arguments (5 given)

comment:63 Changed 3 years ago by jdemeyer

  • Milestone changed from sage-5.1 to sage-5.2

comment:64 Changed 2 years ago by jdemeyer

  • Status changed from positive_review to needs_work

This doesn't apply cleanly to sage-5.1.beta6:

applying /release/merger/patches/trac_11143_v2.5.patch
patching file sage/functions/all.py
Hunk #1 succeeded at 25 with fuzz 1 (offset 0 lines).
patching file sage/interfaces/maxima_lib.py
Hunk #1 FAILED at 1354
1 out of 1 hunks FAILED -- saving rejects to file sage/interfaces/maxima_lib.py.rej
abort: patch failed to apply

comment:65 Changed 2 years ago by kcrisman

This is a tiny thing due to the Lambert W function, #11888. I'll upload a clean one to beta6 momentarily.

Changed 2 years ago by kcrisman

Rebased to 5.1.beta6

comment:66 Changed 2 years ago by kcrisman

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

Patchbot:

Apply to the Sage library:

  1. trac_11143-v2.5-rebased.patch
  2. trac-11143-ref.2.patch

This should be fine.

comment:67 Changed 2 years ago by kcrisman

  • Status changed from positive_review to needs_work

Hmm, I get an error here. This is on 5.1.beta6 with these patches only in addition.

$ ../../sage -t sage/symbolic/pynac.pyx
sage -t  "devel/sage-main/sage/symbolic/pynac.pyx"          
**********************************************************************
File "/Users/.../sage-5.1.beta6/devel/sage-main/sage/symbolic/pynac.pyx", line 544:
    sage: get_sfunction_from_serial(i) == foo
Expected:
    True
Got:
    False
**********************************************************************
File "/Users/.../sage-5.1.beta6/devel/sage-main/sage/symbolic/pynac.pyx", line 546:
    sage: py_latex_fderivative(i, (0, 1, 0, 1), (x, y^z))
Expected:
    D[0, 1, 0, 1]func_with_args(x, y^z)
Got:
    D[0, 1, 0, 1]\left({\rm foo}\right)\left(x, y^{z}\right)
**********************************************************************
1 items had failures:
   2 of  19 in __main__.example_18
***Test Failed*** 2 failures.
For whitespace errors, see the file /Users/.../.sage//tmp/pynac_4902.py
	 [5.7 s]
 
----------------------------------------------------------------------
The following tests failed:


	sage -t  "devel/sage-main/sage/symbolic/pynac.pyx"
Total time for all tests: 5.7 seconds
Last edited 2 years ago by kcrisman (previous) (diff)

comment:68 Changed 2 years ago by kcrisman

I am not sure what is going on here. Testing this by hand works, testing it in the doctest doesn't. It's definitely the main patch that causes this.

I had a few semi-educated guesses, but they didn't seem to work. I'll keep digging.

Last edited 2 years ago by kcrisman (previous) (diff)

comment:69 in reply to: ↑ 52 ; follow-up: Changed 2 years ago by kcrisman

Okay, I've reported the issue with Maxima using the wrong branch for Shi at this bug tracker ticket.

(%i30) float(expintegral_shi(1/2));
(%o30) 3.141592653589793*%i+0.506996749819667

Incidentally, this is apparently now fixed. Is this something we should test for (not in this ticket) as well, since we use mpmath for it? I really hate how there is no link on the Maxima bugtracker as to the actual change in the code. Makes it harder for new folks to jump in.

comment:70 Changed 2 years ago by kcrisman

Aagh! It was what I originally thought - that we just have too many serials now. But I kept on testing the wrong occurrences of this code.

Patch coming up should fix this and at least postpone the agony the next time we add functions...

Changed 2 years ago by kcrisman

comment:71 Changed 2 years ago by kcrisman

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

Patchbot:

Apply to the Sage library:

  1. trac_11143-v2.5-rebased.patch
  2. trac-11143-ref.2.patch
  3. trac_11143-pynac-serials.patch

I feel like the latest patch is from an obscure enough piece of code that it could need review. On the other hand, it is a pretty trivial patch in other ways and only changes doctests. I'll put positive review for now, but feel free to change it back if you feel otherwise, Jeroen; I'm sure burcin or benjaminfjones would be able to sign off on it pretty quickly otherwise.

comment:72 Changed 2 years ago by benjaminfjones

The patch to the doctests in sage/symbolic/pnyac.pyx look (and tests) fine. I suspected it would be something like this, but didn't have the time earlier to track it down. Thanks!

comment:73 in reply to: ↑ 62 Changed 2 years ago by kcrisman

Replying to JoalHeagney:

With the patches applied to 5.0, my fractional integral code works on sine functions, but I still get the following error if pretty printing is on. Is there a trac post for this?

Sorry not to get to this earlier, Joal. I have to admit that I haven't got a clue on this one yet, but don't have time to track it down quickly, since this is your own custom code thus far (though perhaps someday it might be useful for Sage?). Why don't you open a ticket for it, saying something about a pretty print error (as opposed to about fractional integrals). Then at least it has its own place, and you can document exactly what does and doesn't work.

comment:74 in reply to: ↑ 69 ; follow-up: Changed 2 years ago by zimmerma

Replying to kcrisman:

> > (%i30) float(expintegral_shi(1/2));
> > (%o30) 3.141592653589793*%i+0.506996749819667

Is this something we should test for?

yes please!

Paul

comment:75 in reply to: ↑ 74 Changed 2 years ago by kcrisman

> > (%i30) float(expintegral_shi(1/2));
> > (%o30) 3.141592653589793*%i+0.506996749819667

Is this something we should test for?

yes please!

Keep in mind this is a Maxima output, not one we currently can even easily access without using Maxima, since this ticket uses mpmath to evaluate these functions, and we don't (yet) have different algorithm parameters for all these special functions. We do test a fair number of values of the "Sage" versions of these functions in this ticket!

comment:76 Changed 2 years ago by vbraun

  • Dependencies set to #13109
  • Status changed from positive_review to needs_work

Switching the deprecation to the new syntax.

comment:77 Changed 2 years ago by vbraun

  • Authors changed from Benjamin Jones to Benjamin Jones, Volker Braun
  • Description modified (diff)
  • Status changed from needs_work to positive_review

The new patch is just rediffed to apply on top of #13109. Switching back to positive review.

comment:78 Changed 2 years ago by jdemeyer

  • Milestone changed from sage-5.2 to sage-pending

Changed 2 years ago by vbraun

Rediffed patch

comment:79 Changed 2 years ago by vbraun

The authors of #11475 thought its a good idea to remove trailing whitespace in random places, which broke this patch. Rediffed for sage-5.2.beta1.

comment:80 Changed 2 years ago by kcrisman

What is this pending? If #11475 has been unmerged, then is there an appropriate patch that can be merged in Sage 5.3, since 5.2.rc0 is out?

comment:81 Changed 2 years ago by kcrisman

  • Description modified (diff)

AND the patchbot says it doesn't work. So there. Let's see if this works:

Apply to the Sage library:

  1. trac_11143-v2.5-rebased.patch
  2. trac-11143-ref.2.patch
  3. trac_11143-pynac-serials.patch

comment:82 Changed 2 years ago by kcrisman

  • Milestone changed from sage-pending to sage-5.3

All tests in sage/calculus, sage/symbolic, and sage/functions pass against 5.2.beta1. If patchbot later says it's ok, then I think this is ok. Puttting to Sage 5.3.

comment:83 Changed 2 years ago by jdemeyer

It was put to sage-pending because of #13109 (which is now merged).

comment:84 follow-up: Changed 2 years ago by kcrisman

Oh great. Unfortunately, Volker apparently replaced the patch based on #13109 with one based on #11475, and the previous rebase was before #13109, so I don't think there is any patch that would currently apply... this is silliness. I thought that #13109 had been delayed because of the sage-combinat issue, but apparently not?

comment:85 Changed 2 years ago by kcrisman

Hey Ben, question. Which Li/li is Maxima using with its conversion? I think that expintegral_li actually correspond to Li (#3401), not the li provided here, though the integral of log_integral does seem right... Also, while looking at this I noticed that

\operatorname{li}(x) = \int_0^z \frac{dt}{\ln(t)} = \operatorname{Ei}(\ln(x))

probably should have a x not a z in the integral bound.

comment:86 in reply to: ↑ 84 Changed 2 years ago by jdemeyer

  • Status changed from positive_review to needs_work

Replying to kcrisman:

Oh great. Unfortunately, Volker apparently replaced the patch based on #13109 with one based on #11475, and the previous rebase was before #13109, so I don't think there is any patch that would currently apply...

Indeed, this needs to be rebased to sage-5.2.rc1.

Changed 2 years ago by vbraun

Rebased patch

comment:87 Changed 2 years ago by vbraun

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

The rebased patch applies on sage-5.2.rc1. Switching this back to positive review

comment:88 Changed 2 years ago by jdemeyer

  • Status changed from positive_review to needs_work

There is a trivial doctest failure on 32-bit i386 systems (possibly other systems too):

sage -t  --long -force_lib devel/sage/sage/functions/exp_integral.py
**********************************************************************
File "/var/lib/buildbot/build/sage/arando-1/arando_full/build/sage-5.3.beta1/devel/sage-main/sage/functions/exp_integral.py", line 1297:
    sage: w = exponential_integral_1(2,4); w
Expected:
    [0.04890051070806112, 0.0037793524098489067, 0.00036008245216265873, 3.7665622843924751e-05] 
Got:
    [0.04890051070806112, 0.0037793524098489067, 0.00036008245216265873, 3.766562284392475e-05]
**********************************************************************

Changed 2 years ago by benjaminfjones

fixed failing doctest

comment:89 Changed 2 years ago by benjaminfjones

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

Fixed the failing doctest (for the reviewer: this only changed lines 1300--1301 in sage/functions/exp_integral.py).

comment:90 Changed 2 years ago by jdemeyer

  • Status changed from needs_review to needs_work

Wouldn't a relative tolerance be better than an absolute one? With the absolute tolerance, the fourth value can essentially be anything.

comment:91 Changed 2 years ago by benjaminfjones

Yes, you're right. I'll change that.

Changed 2 years ago by benjaminfjones

changed abs tol to rel tol

comment:92 Changed 2 years ago by benjaminfjones

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

comment:93 Changed 2 years ago by benjaminfjones

Patchbot, please kindly apply:

to the Sage library.

Thank you in advance.

comment:94 Changed 2 years ago by jdemeyer

I'm happy with the changes to that line. I haven't re-checked the whole patch again, but I thrust that you didn't change anything else.

comment:95 Changed 2 years ago by jdemeyer

  • Status changed from needs_review to positive_review

comment:96 Changed 2 years ago by jdemeyer

  • Merged in set to sage-5.3.beta2
  • Resolution set to fixed
  • Status changed from positive_review to closed

comment:97 Changed 18 months ago by kcrisman

See #14766 for a followup.

Note: See TracTickets for help on using tickets.