Opened 12 years ago
Last modified 9 years ago
#9706 closed enhancement
New Version of orthogonal Polynomials — at Version 28
Reported by: | maldun | Owned by: | burcin, maldun |
---|---|---|---|
Priority: | major | Milestone: | sage-6.1 |
Component: | symbolics | Keywords: | orthogonal polynomials, symbolics |
Cc: | fredrik.johansson, fstan, kcrisman | Merged in: | |
Authors: | Stefan Reiterer | Reviewers: | |
Report Upstream: | N/A | Work issues: | |
Branch: | Commit: | ||
Dependencies: | Stopgaps: |
Description (last modified by )
The current implementation of orthogonal polynomials is just a wrapper around maxima. (see http://wiki.sagemath.org/symbolics/) This update holds the following changes:
*using of the pynac class for symbolic functions. *faster evaluation in general *evaluation of special values *mpmath for numeric evaluation
Remarks:
*The current patch needs scipy-0.8. One has to install it before testing (see #9808 for spkg's and installation instructions)
- Some of the old doctests in the old file don't work any more, due to coercion problems with pynac (see #9769)
- Some doctests in Sage change, due to the fact that new BuiltIn? functions are added. symbolic.random_test.py had output changes since the random expression creation changed of course. The tests in pynac.pyx also changed, but this has a strange behavior (see below).
Change History (38)
Changed 12 years ago by
Changed 12 years ago by
Newer version, with legendre_P, and faster evaluation of symbolic expressions
comment:1 follow-up: ↓ 2 Changed 12 years ago by
All Polys now have their own class. Much faster evaluation is added. Numerical evaluation is provided. Except for legendre_Q, gen_legendre_P, and gen_legendre_Q these aren't ready yet
comment:2 in reply to: ↑ 1 Changed 12 years ago by
Replying to maldun:
All Polys now have their own class. Much faster evaluation is added. Numerical evaluation is provided. Except for legendre_Q, gen_legendre_P, and gen_legendre_Q these aren't ready yet
orthogonal_polys4.py hold all changes but is not a patch yet, because it holds old code fragments, which I have to clean up...
comment:3 Changed 12 years ago by
- Cc fredrik.johansson added
comment:4 Changed 12 years ago by
I added in the latest patch (and orthogonal_polys.4.py contains these changes also) a new symbolic evaluation method for the orthogonal polynomials: Instead of call Maxima or use of the recursion, the polynomial is evaluated just using explicit formulas from Abramowitz and Stegun. This is an O(n) algorithm of course.
a little comparison on my machine: old version:
sage: time chebyshev_T(10,x); CPU times: user 0.01 s, sys: 0.00 s, total: 0.01 s Wall time: 0.04 s sage: time chebyshev_T(100,x); CPU times: user 0.13 s, sys: 0.01 s, total: 0.14 s Wall time: 0.23 s sage: time chebyshev_T(1000,x); CPU times: user 5.01 s, sys: 0.01 s, total: 5.02 s Wall time: 6.98 s sage time chebyshev_T(5000,x); ??? (I got no output her after 2min)
sage: time gegenbauer(10,5,x); CPU times: user 0.01 s, sys: 0.00 s, total: 0.01 s Wall time: 0.05 s sage: time gegenbauer(100,5,x); CPU times: user 0.19 s, sys: 0.00 s, total: 0.19 s Wall time: 0.29 s sage: time gegenbauer(1000,5,x); CPU times: user 5.46 s, sys: 0.02 s, total: 5.48 s Wall time: 7.79 s
New Version sage: time chebyshev_T(10,x); CPU times: user 0.01 s, sys: 0.00 s, total: 0.01 s Wall time: 0.01 s sage: time chebyshev_T(100,x); CPU times: user 0.06 s, sys: 0.00 s, total: 0.06 s Wall time: 0.08 s sage: time chebyshev_T(1000,x); CPU times: user 1.22 s, sys: 0.00 s, total: 1.22 s Wall time: 1.22 s sage: time chebyshev_T(5000,x); CPU times: user 27.17 s, sys: 0.15 s, total: 27.32 s Wall time: 27.46 s
sage: time gegenbauer(10,5,x); CPU times: user 0.01 s, sys: 0.00 s, total: 0.01 s Wall time: 0.01 s sage: time gegenbauer(100,5,x); CPU times: user 0.03 s, sys: 0.00 s, total: 0.03 s Wall time: 0.04 s sage: time gegenbauer(1000,5,x); CPU times: user 1.08 s, sys: 0.01 s, total: 1.09 s Wall time: 1.11 s
A little bit faster :) I also don't need to spawn an instance of maxima which makes the initialisation faster.
And now also wider symbolic evaluation is possible:
old version: sage: var('a') a sage: gegenbauer(3,a,x) ... NameError?: name 'a' is not defined
new version: sage: var('a') a sage: gegenbauer(3,a,x) 4/3*x3*gamma(a + 3) - 2*x*gamma(a + 2)
The code needs now some cleanup, especially the documentations. The complete versions for legendre_Q, gen_legendre_P, and gen_legendre_Q will not be finished soon since the mpmath functions, don't seem to work correctly... I only provide a call function for maxima for them now.
comment:5 follow-up: ↓ 7 Changed 12 years ago by
The complete versions for legendre_Q, gen_legendre_P, and gen_legendre_Q will not be finished soon since the mpmath functions, don't seem to work correctly...
Care to elaborate?
comment:6 Changed 12 years ago by
Killed bug in legendre_P
comment:7 in reply to: ↑ 5 Changed 12 years ago by
Replying to fredrik.johansson:
The complete versions for legendre_Q, gen_legendre_P, and gen_legendre_Q will not be finished soon since the mpmath functions, don't seem to work correctly...
Care to elaborate?
Sorry for the late answer, I was on holidays.
In mpmath I have probs with the legenp and legenq functions. For some inputs I get this error:
sage: mpmath.call(mpmath.legenp,5,1,2) --------------------------------------------------------------------------- OverflowError Traceback (most recent call last) /home/maldun/prog/sage/ortho/<ipython console> in <module>() /home/maldun/sage/sage-4.5.1/local/lib/python2.6/site-packages/sage/libs/mpmath/utils.so in sage.libs.mpmath.utils.call (sage/libs/mpmath/utils.c:5021)() /home/maldun/sage/sage-4.5.1/local/lib/python2.6/site-packages/mpmath/functions/hypergeometric.pyc in legenp(ctx, n, m, z, type, **kwargs) 1481 T = [1+z, 1-z], [g, -g], [], [1-m], [-n, n+1], [1-m], 0.5*(1-z) 1482 return (T,) -> 1483 return ctx.hypercomb(h, [n,m], **kwargs) 1484 if type == 3: 1485 def h(n,m): /home/maldun/sage/sage-4.5.1/local/lib/python2.6/site-packages/mpmath/functions/hypergeometric.pyc in hypercomb(ctx, function, params, discard_known_zeros, **kwargs) 125 [ctx.gamma(a) for a in alpha_s] + \ 126 [ctx.rgamma(b) for b in beta_s] + \ --> 127 [ctx.power(w,c) for (w,c) in zip(w_s,c_s)]) 128 if verbose: 129 print " Value:", v /home/maldun/sage/sage-4.5.1/local/lib/python2.6/site-packages/mpmath/ctx_base.pyc in power(ctx, x, y) 417 3.16470269330255923143453723949e+12978188 418 """ --> 419 return ctx.convert(x) ** ctx.convert(y) 420 421 def _zeta_int(ctx, n): /home/maldun/sage/sage-4.5.1/local/lib/python2.6/site-packages/sage/libs/mpmath/ext_main.so in sage.libs.mpmath.ext_main.mpnumber.__pow__ (sage/libs/mpmath/ext_main.c:13946)() /home/maldun/sage/sage-4.5.1/local/lib/python2.6/site-packages/sage/libs/mpmath/ext_main.so in sage.libs.mpmath.ext_main.binop (sage/libs/mpmath/ext_main.c:4588)() /home/maldun/sage/sage-4.5.1/local/lib/python2.6/site-packages/mpmath/libmp/libelefun.pyc in mpf_pow(s, t, prec, rnd) 340 # General formula: s**t = exp(t*log(s)) 341 # TODO: handle rnd direction of the logarithm carefully --> 342 c = mpf_log(s, prec+10, rnd) 343 return mpf_exp(mpf_mul(t, c), prec, rnd) 344 /home/maldun/sage/sage-4.5.1/local/lib/python2.6/site-packages/mpmath/libmp/libelefun.pyc in mpf_log(x, prec, rnd) 725 # optimal between 1000 and 100,000 digits. 726 if wp <= LOG_TAYLOR_PREC: --> 727 m = log_taylor_cached(lshift(man, wp-bc), wp) 728 if mag: 729 m += mag*ln2_fixed(wp) /home/maldun/sage/sage-4.5.1/local/lib/python2.6/site-packages/mpmath/libmp/libelefun.pyc in log_taylor_cached(x, prec) 643 else: 644 a = n << (cached_prec - LOG_TAYLOR_SHIFT) --> 645 log_a = log_taylor(a, cached_prec, 8) 646 log_taylor_cache[n, cached_prec] = (a, log_a) 647 a >>= dprec /home/maldun/sage/sage-4.5.1/local/lib/python2.6/site-packages/mpmath/libmp/libelefun.pyc in log_taylor(x, prec, r) 607 """ 608 for i in xrange(r): --> 609 x = isqrt_fast(x<<prec) 610 one = MPZ_ONE << prec 611 v = ((x-one)<<prec)//(x+one) /home/maldun/sage/sage-4.5.1/local/lib/python2.6/site-packages/mpmath/libmp/libintmath.pyc in isqrt_fast_python(x) 240 y = (y + x//y) >> 1 241 return y --> 242 bc = bitcount(x) 243 guard_bits = 10 244 x <<= 2*guard_bits /home/maldun/sage/sage-4.5.1/local/lib/python2.6/site-packages/mpmath/libmp/libintmath.pyc in python_bitcount(n) 78 if bc != 300: 79 return bc ---> 80 bc = int(math.log(n, 2)) - 4 81 return bc + bctable[n>>bc] 82 OverflowError: cannot convert float infinity to integer
comment:8 follow-ups: ↓ 9 ↓ 10 Changed 12 years ago by
That looks strange. I get:
sage: import sage.libs.mpmath.all as mpmath sage: mpmath.call(mpmath.legenp, 5,1,2) -2.96434298694874e-22 - 912.574269237852*I sage: mpmath.call(mpmath.legenp, 5,1,2, prec=100) -2.1062923756778274648015607872e-36 - 912.57426923785222402727329118*I
comment:9 in reply to: ↑ 8 Changed 12 years ago by
Replying to fredrik.johansson:
That looks strange. I get:
sage: import sage.libs.mpmath.all as mpmath sage: mpmath.call(mpmath.legenp, 5,1,2) -2.96434298694874e-22 - 912.574269237852*I sage: mpmath.call(mpmath.legenp, 5,1,2, prec=100) -2.1062923756778274648015607872e-36 - 912.57426923785222402727329118*I
Hm strange. Today I install the new Sage version, perhaps it will then work again
comment:10 in reply to: ↑ 8 Changed 12 years ago by
Replying to fredrik.johansson:
That looks strange. I get:
sage: import sage.libs.mpmath.all as mpmath sage: mpmath.call(mpmath.legenp, 5,1,2) -2.96434298694874e-22 - 912.574269237852*I sage: mpmath.call(mpmath.legenp, 5,1,2, prec=100) -2.1062923756778274648015607872e-36 - 912.57426923785222402727329118*I
It was the old version!a Thanx for pointing that out, I will continue soon =)
comment:11 follow-ups: ↓ 12 ↓ 14 Changed 12 years ago by
So now a "beta" is ready with full support of all classes.
Only the Legendre functions are still using Maxima.
some advances for the future:
-Zernike polys (this should be done in the next time, since explicit formulas are available) -support for numpy_eval. (But this will be done, when the scipy package is updated to 0.8, else it has no sense, because the current version of scipy does not support ortho polys well, but the newer can handle them)
Now I need some people for testing this out =)
comment:12 in reply to: ↑ 11 Changed 12 years ago by
And there was an interisting bug:
the import of mpmath at the beginning of the file caused the whole trouble I had with the numeric evaluation of the legendre functions....
I think I should report this..
comment:13 Changed 12 years ago by
- Type changed from defect to enhancement
comment:14 in reply to: ↑ 11 Changed 12 years ago by
-support for numpy_eval. (But this will be done, when the scipy package is updated to 0.8, else it has no sense, because the current version of scipy does not support ortho polys well, but the newer can handle them)
I decided to give at least some numpy support for compability reasons. But this is a bad hack...when scipy 0.8 comes I use scipy itself, I change this to a better version :)
comment:15 Changed 12 years ago by
- Status changed from new to needs_review
comment:16 Changed 12 years ago by
- Milestone set to sage-5.0
comment:17 Changed 12 years ago by
Some of the old doctests fail. But it is not my fault, it seem's that it is a bug in the SymbolicFunction? class.
comment:18 Changed 12 years ago by
- Milestone changed from sage-5.0 to sage-4.5.3
comment:19 Changed 12 years ago by
- Owner changed from burcin to burcin, maldun
comment:20 follow-up: ↓ 21 Changed 12 years ago by
Hi Stefan,
can you post a patch corresponding to attachment:orthogonal_polys.8.py for review?
Thanks,
Burcin
comment:21 in reply to: ↑ 20 Changed 12 years ago by
Replying to burcin:
Hi Stefan,
can you post a patch corresponding to attachment:orthogonal_polys.8.py for review?
Thanks,
Burcin
Done!
comment:22 follow-up: ↓ 23 Changed 12 years ago by
Why is mpmath's precision used by default? Shouldn't the default be RR / CC precision? Actually, does _evalf_ ever get called without this information?
Some complex tests would be nice.
comment:23 in reply to: ↑ 22 ; follow-up: ↓ 24 Changed 12 years ago by
Replying to fredrik.johansson:
Why is mpmath's precision used by default? Shouldn't the default be RR / CC precision? Actually, does _evalf_ ever get called without this information?
Some complex tests would be nice.
This is a good point, and it shouldn't be a problem to change that. But I don't think it's a big deal, because the function takes the "parents" precision, which means, if my input is RR it evals it with RR's precision.
Of course can you call _evalf_ just with (), and then the default value is used.
I just sticked to the old's version tests, and expanded it. Of course it's possible to expand the tests. I hope I will find some time for it soon, since I have some other more urgent things todo also.
comment:24 in reply to: ↑ 23 Changed 12 years ago by
Replying to maldun:
Replying to fredrik.johansson:
Why is mpmath's precision used by default? Shouldn't the default be RR / CC precision? Actually, does _evalf_ ever get called without this information?
Some complex tests would be nice.
This is a good point, and it shouldn't be a problem to change that. But I don't think it's a big deal, because the function takes the "parents" precision, which means, if my input is RR it evals it with RR's precision.
Of course can you call _evalf_ just with (), and then the default value is used.
Ok sorry, wrong explination: when your input are exact data types like ZZ ore QQ then the parent has no precision, then you need a default value
comment:25 Changed 12 years ago by
Since it seems that numpy-1.4.1, and scipy 0.8 should work now (see #9808) I programmed a version which uses scipy itself to evaluate the orthogonal polys for numpy arrays. When the new versions of numpy/scipy become merged into sage I will provide a patch for these.
Another thing I have to mention are these 2 failde doctests:
- sage -t -long "devel/sage/sage/symbolic/random_tests.py"
- sage -t -long "devel/sage/sage/symbolic/pynac.pyx"
sage -t -long "devel/sage/sage/symbolic/random_tests.py" ********************************************************************** File "/home/maldun/sage/sage-4.5.2/devel/sage/sage/symbolic/random_tests.py", line 17: sage: [f for (one,f,arity) in _mk_full_functions()] Expected: [Ei, abs, arccos, arccosh, arccot, arccoth, arccsc, arccsch, arcsec, arcsech, arcsin, arcsinh, arctan, arctan2, arctanh, binomial, ceil, conjugate, cos, cosh, cot, coth, csc, csch, dickman_rho, dilog, dirac_delta, elliptic_e, elliptic_ec, elliptic_eu, elliptic_f, elliptic_kc, elliptic_pi, erf, exp, factorial, floor, heaviside, imag_part, integrate, kronecker_delta, log, polylog, real_part, sec, sech, sgn, sin, sinh, tan, tanh, unit_step, zeta, zetaderiv] Got: [Ei, abs, arccos, arccosh, arccot, arccoth, arccsc, arccsch, arcsec, arcsech, arcsin, arcsinh, arctan, arctan2, arctanh, binomial, ceil, chebyshev_T, chebyshev_U, conjugate, cos, cosh, cot, coth, csc, csch, dickman_rho, dilog, dirac_delta, elliptic_e, elliptic_ec, elliptic_eu, elliptic_f, elliptic_kc, elliptic_pi, erf, exp, factorial, floor, gegenbauer, gen_laguerre, gen_legendre_P, gen_legendre_Q, heaviside, hermite, imag_part, integrate, jacobi_P, kronecker_delta, laguerre, legendre_P, legendre_Q, log, polylog, real_part, sec, sech, sgn, sin, sinh, tan, tanh, unit_step, zeta, zetaderiv] ********************************************************************** File "/home/maldun/sage/sage-4.5.2/devel/sage/sage/symbolic/random_tests.py", line 237: sage: random_expr(50, nvars=3, coeff_generator=CDF.random_element) Expected: (euler_gamma - v3^(-e) + (v2 - factorial(-e/v2))^(((2.85879036573 - 1.18163393202*I)*v2 + (2.85879036573 - 1.18163393202*I)*v3)*pi - 0.247786879678 + 0.931826724898*I)*arccsc((0.891138386848 - 0.0936820840629*I)/v1) + (-0.553423153995 + 0.5481180572*I)*v3 + 0.149683576515 - 0.155746451854*I)*v1 + arccsch(pi + e)*elliptic_f(khinchin*v2, 1.4656989704 + 0.863754357069*I) Got: -v1*e^((0.0666829501658 + 0.206976992303*I)/(v3 + e))/v3 + hermite(-(v3^(-0.48519994364 - 0.485764091302*I) - log((1.21734510331 - 1.22580558833*I)*pi*v1 + zeta((0.781366128261 + 0.957400336147*I)*v1*e + (-1.8919687109 + 0.753422167447*I)*elliptic_f(v1, v1))*arccsch(v3)))*v1, (-0.647983235144 + 1.20665952957*I)*v1 + (0.0909404921682 + 0.281538203756*I)/v3) ********************************************************************** File "/home/maldun/sage/sage-4.5.2/devel/sage/sage/symbolic/random_tests.py", line 239: sage: random_expr(5, verbose=True) Exception raised: Traceback (most recent call last): File "/home/maldun/sage/sage-4.5.2/local/bin/ncadoctest.py", line 1231, in run_one_test self.run_one_example(test, example, filename, compileflags) File "/home/maldun/sage/sage-4.5.2/local/bin/sagedoctest.py", line 38, in run_one_example OrigDocTestRunner.run_one_example(self, test, example, filename, compileflags) File "/home/maldun/sage/sage-4.5.2/local/bin/ncadoctest.py", line 1172, in run_one_example compileflags, 1) in test.globs File "<doctest __main__.example_5[5]>", line 1, in <module> random_expr(Integer(5), verbose=True)###line 239: sage: random_expr(5, verbose=True) File "/home/maldun/sage/sage-4.5.2/local/lib/python/site-packages/sage/symbolic/random_tests.py", line 254, in random_expr return random_expr_helper(size, internal, leaves, verbose) File "/home/maldun/sage/sage-4.5.2/local/lib/python/site-packages/sage/symbolic/random_tests.py", line 210, in random_expr_helper return r[1](*children) File "element.pyx", line 1529, in sage.structure.element.RingElement.__div__ (sage/structure/element.c:11992) File "coerce.pyx", line 713, in sage.structure.coerce.CoercionModel_cache_maps.bin_op (sage/structure/coerce.c:6126) File "element.pyx", line 1527, in sage.structure.element.RingElement.__div__ (sage/structure/element.c:11973) File "expression.pyx", line 2269, in sage.symbolic.expression.Expression._div_ (sage/symbolic/expression.cpp:11444) ZeroDivisionError: Symbolic division by zero ********************************************************************** 2 items had failures: 1 of 4 in __main__.example_0 2 of 6 in __main__.example_5 ***Test Failed*** 3 failures. For whitespace errors, see the file /home/maldun/.sage//tmp/.doctest_random_tests.py [7.7 s] ---------------------------------------------------------------------- The following tests failed: sage -t -long "devel/sage/sage/symbolic/random_tests.py" Total time for all tests: 7.8 seconds
I quite understand these, because we have introduced new functions, but I don't understand the exception in the last one
sage -t -long "devel/sage/sage/symbolic/pynac.pyx" ********************************************************************** File "/home/maldun/sage/sage-4.5.2/devel/sage/sage/symbolic/pynac.pyx", line 386: sage: get_sfunction_from_serial(i) == foo Expected: True Got: False ********************************************************************** File "/home/maldun/sage/sage-4.5.2/devel/sage/sage/symbolic/pynac.pyx", line 388: sage: py_latex_function_pystring(i, (x,y^z)) Expected: 'my args are: x, y^z' Got: '\\mathrm{bar}\\left(x, y^{z}\\right)' ********************************************************************** File "/home/maldun/sage/sage-4.5.2/devel/sage/sage/symbolic/pynac.pyx", line 478: sage: get_sfunction_from_serial(i) == foo Expected: True Got: False ********************************************************************** File "/home/maldun/sage/sage-4.5.2/devel/sage/sage/symbolic/pynac.pyx", line 480: sage: py_print_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](foo)(x, y^z) ********************************************************************** File "/home/maldun/sage/sage-4.5.2/devel/sage/sage/symbolic/pynac.pyx", line 540: sage: get_sfunction_from_serial(i) == foo Expected: True Got: False ********************************************************************** File "/home/maldun/sage/sage-4.5.2/devel/sage/sage/symbolic/pynac.pyx", line 542: 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(\mathrm{bar}\right)\left(x, y^{z}\right) ********************************************************************** 3 items had failures: 2 of 19 in __main__.example_14 2 of 14 in __main__.example_16 2 of 18 in __main__.example_18 ***Test Failed*** 6 failures. For whitespace errors, see the file /home/maldun/.sage//tmp/.doctest_pynac.py [7.3 s] ---------------------------------------------------------------------- The following tests failed: sage -t -long "devel/sage/sage/symbolic/pynac.pyx" Total time for all tests: 7.3 seconds
And these are really strange, because when I type then into sage by hand everything works. wtf?? Can anyone have a look at these?
comment:26 Changed 12 years ago by
- Milestone changed from sage-4.6 to sage-5.0
- Status changed from needs_review to needs_work
comment:27 Changed 12 years ago by
Just cc:ing myself by commenting.
Also, there seems to be a lot of stuff in the latest Python file that is the same as the original one (in terms of explanation, not code). Maybe posting an updated patch (once the numpy/scipy-fest is over, which is hopefully the case) would help some of us figure this out. Thanks for working on this - there is still a lot of overhauling that symbolics could use, but this is a great step.
comment:28 Changed 12 years ago by
- Description modified (diff)
- Status changed from needs_work to needs_review
@kcrisman thanks for paying attention. I added now an updated patch and extended instructions.
the doctest changes in symbolic.random_tests.py are easy to explain: new functions are involved -> new random expressions. But I had to change
random_expr(50, nvars=3, coeff_generator=CDF.random_element)
to random_expr(60, nvars=3, coeff_generator=CDF.random_element)
or else one gets an expression generated where a division through zero occours.
As mentioned on sage-devel I repaired the doctests in symbolic.pynac.pyx, the trick is to enlarge the range of the for loop:
for i in range(get_ginac_serial(), get_ginac_serial()+50):
changed to
for i in range(get_ginac_serial(), get_ginac_serial()+100):
now it works. My explaination: since we have new functions we have longer to search, and then we reach our goal. What I can not explain is, that it works, when I type it in by hand.
All doctests pass now, so I think a review would be nice.
-maldun
A new version of the orthogonal_polys.py file.