FAIL: test_yn_zeros, test_ynp_zeros, Negative-order Bessels
Hi all, I do get three failures wrt. to special functions and python2.6 ====================================================================== FAIL: test_yn_zeros (test_basic.TestBessel) ---------------------------------------------------------------------- Traceback (most recent call last): File "/home/nwagner/local/lib64/python2.6/site-packages/scipy/special/tests/test_basic.py", line 1598, in test_yn_zeros 488.98055964441374646], rtol=1e-19) File "/home/nwagner/local/lib64/python2.6/site-packages/scipy/special/tests/test_basic.py", line 38, in assert_tol_equal verbose=verbose, header=header) File "/home/nwagner/local/lib64/python2.6/site-packages/numpy/testing/utils.py", line 295, in assert_array_compare raise AssertionError(msg) AssertionError: Not equal to tolerance rtol=1e-19, atol=0 (mismatch 100.0%) x: array([ 450.136, 463.057, 472.807, 481.274, 488.981]) y: array([ 450.136, 463.057, 472.807, 481.274, 488.981]) ====================================================================== FAIL: test_ynp_zeros (test_basic.TestBessel) ---------------------------------------------------------------------- Traceback (most recent call last): File "/home/nwagner/local/lib64/python2.6/site-packages/scipy/special/tests/test_basic.py", line 1604, in test_ynp_zeros assert_tol_equal(yvp(443, ao), 0, atol=1e-15) File "/home/nwagner/local/lib64/python2.6/site-packages/scipy/special/tests/test_basic.py", line 38, in assert_tol_equal verbose=verbose, header=header) File "/home/nwagner/local/lib64/python2.6/site-packages/numpy/testing/utils.py", line 295, in assert_array_compare raise AssertionError(msg) AssertionError: Not equal to tolerance rtol=1e-07, atol=1e-15 (mismatch 100.0%) x: array([ 1.239e-10, -8.119e-16, 3.608e-16, 5.898e-16, 1.226e-15]) y: array(0) ====================================================================== FAIL: Negative-order Bessels ---------------------------------------------------------------------- Traceback (most recent call last): File "/home/nwagner/local/lib64/python2.6/site-packages/nose-0.10.4-py2.6.egg/nose/case.py", line 182, in runTest self.test(*self.arg) File "/home/nwagner/local/lib64/python2.6/site-packages/scipy/special/tests/test_basic.py", line 1665, in test_ticket_853 assert_tol_equal(iv(-0.5, 1 ), 1.231200214592967) File "/home/nwagner/local/lib64/python2.6/site-packages/scipy/special/tests/test_basic.py", line 38, in assert_tol_equal verbose=verbose, header=header) File "/home/nwagner/local/lib64/python2.6/site-packages/numpy/testing/utils.py", line 295, in assert_array_compare raise AssertionError(msg) AssertionError: Not equal to tolerance rtol=1e-07, atol=0 (mismatch 100.0%) x: array(0.0) y: array(1.231200214592967) Nils
On Sun, 25 Jan 2009 20:03:26 +0100 "Nils Wagner" <nwagner@iam.uni-stuttgart.de> wrote:
Hi all,
I do get three failures wrt. to special functions and python2.6
====================================================================== FAIL: test_yn_zeros (test_basic.TestBessel) ---------------------------------------------------------------------- Traceback (most recent call last): File "/home/nwagner/local/lib64/python2.6/site-packages/scipy/special/tests/test_basic.py", line 1598, in test_yn_zeros 488.98055964441374646], rtol=1e-19) File "/home/nwagner/local/lib64/python2.6/site-packages/scipy/special/tests/test_basic.py", line 38, in assert_tol_equal verbose=verbose, header=header) File "/home/nwagner/local/lib64/python2.6/site-packages/numpy/testing/utils.py", line 295, in assert_array_compare raise AssertionError(msg) AssertionError: Not equal to tolerance rtol=1e-19, atol=0
(mismatch 100.0%) x: array([ 450.136, 463.057, 472.807, 481.274, 488.981]) y: array([ 450.136, 463.057, 472.807, 481.274, 488.981])
====================================================================== FAIL: test_ynp_zeros (test_basic.TestBessel) ---------------------------------------------------------------------- Traceback (most recent call last): File "/home/nwagner/local/lib64/python2.6/site-packages/scipy/special/tests/test_basic.py", line 1604, in test_ynp_zeros assert_tol_equal(yvp(443, ao), 0, atol=1e-15) File "/home/nwagner/local/lib64/python2.6/site-packages/scipy/special/tests/test_basic.py", line 38, in assert_tol_equal verbose=verbose, header=header) File "/home/nwagner/local/lib64/python2.6/site-packages/numpy/testing/utils.py", line 295, in assert_array_compare raise AssertionError(msg) AssertionError: Not equal to tolerance rtol=1e-07, atol=1e-15
(mismatch 100.0%) x: array([ 1.239e-10, -8.119e-16, 3.608e-16, 5.898e-16, 1.226e-15]) y: array(0)
====================================================================== FAIL: Negative-order Bessels ---------------------------------------------------------------------- Traceback (most recent call last): File "/home/nwagner/local/lib64/python2.6/site-packages/nose-0.10.4-py2.6.egg/nose/case.py", line 182, in runTest self.test(*self.arg) File "/home/nwagner/local/lib64/python2.6/site-packages/scipy/special/tests/test_basic.py", line 1665, in test_ticket_853 assert_tol_equal(iv(-0.5, 1 ), 1.231200214592967) File "/home/nwagner/local/lib64/python2.6/site-packages/scipy/special/tests/test_basic.py", line 38, in assert_tol_equal verbose=verbose, header=header) File "/home/nwagner/local/lib64/python2.6/site-packages/numpy/testing/utils.py", line 295, in assert_array_compare raise AssertionError(msg) AssertionError: Not equal to tolerance rtol=1e-07, atol=0
(mismatch 100.0%) x: array(0.0) y: array(1.231200214592967)
Nils _______________________________________________ Scipy-dev mailing list Scipy-dev@scipy.org http://projects.scipy.org/mailman/listinfo/scipy-dev
scipy.__version__ '0.8.0.dev5524'
====================================================================== FAIL: test_yv_cephes_vs_amos (test_basic.TestBessel) ---------------------------------------------------------------------- Traceback (most recent call last): File "/home/nwagner/local/lib64/python2.6/site-packages/scipy/special/tests/test_basic.py", line 1650, in test_yv_cephes_vs_amos self.check_cephes_vs_amos(yv, yn, rtol=1e-11, atol=1e-305) File "/home/nwagner/local/lib64/python2.6/site-packages/scipy/special/tests/test_basic.py", line 1640, in check_cephes_vs_amos assert c2.imag != 0, (v, z) AssertionError: (301, 1.0) Nils
Sun, 25 Jan 2009 20:25:16 +0100, Nils Wagner wrote: [clip]
scipy.__version__ '0.8.0.dev5524' ====================================================================== FAIL: test_yv_cephes_vs_amos (test_basic.TestBessel)
[clip] Traceback (most recent call last): File "/home/nwagner/local/lib64/python2.6/site-packages/scipy/special/tests/ test_basic.py", line 1650, in test_yv_cephes_vs_amos self.check_cephes_vs_amos(yv, yn, rtol=1e-11, atol=1e-305) File "/home/nwagner/local/lib64/python2.6/site-packages/scipy/special/tests/ test_basic.py", line 1640, in check_cephes_vs_amos assert c2.imag != 0, (v, z) AssertionError: (301, 1.0)
This should have been fixed in r5524. Note that you need to touch scipy/special/_cephesmodule.c before rebuilding, to make sure that scipy.special._cephes module is rebuilt. The build system apparently does not track dependencies properly here. -- Pauli Virtanen
On Mon, 26 Jan 2009 09:31:31 +0000 (UTC) Pauli Virtanen <pav@iki.fi> wrote:
Sun, 25 Jan 2009 20:25:16 +0100, Nils Wagner wrote: [clip]
scipy.__version__ '0.8.0.dev5524' ====================================================================== FAIL: test_yv_cephes_vs_amos (test_basic.TestBessel)
[clip] Traceback (most recent call last): File "/home/nwagner/local/lib64/python2.6/site-packages/scipy/special/tests/ test_basic.py", line 1650, in test_yv_cephes_vs_amos self.check_cephes_vs_amos(yv, yn, rtol=1e-11, atol=1e-305) File "/home/nwagner/local/lib64/python2.6/site-packages/scipy/special/tests/ test_basic.py", line 1640, in check_cephes_vs_amos assert c2.imag != 0, (v, z) AssertionError: (301, 1.0)
This should have been fixed in r5524. Note that you need to
touch scipy/special/_cephesmodule.c
Done, but the failure persists.
before rebuilding, to make sure that scipy.special._cephes module is rebuilt. The build system apparently does not track dependencies properly here.
-- Pauli Virtanen
I am using 64-bit linux CentOS release 4.6 (Final). Python 2.5.1 (r251:54863, Dec 21 2007, 09:21:07) [GCC 3.4.6 20060404 (Red Hat 3.4.6-3)] on linux2 Type "help", "copyright", "credits" or "license" for more information.
import scipy.special scipy.special.test() Running unit tests for scipy.special NumPy version 1.3.0.dev6333 NumPy is installed in /data/home/nwagner/local/lib/python2.5/site-packages/numpy SciPy version 0.8.0.dev5525 SciPy is installed in /data/home/nwagner/local/lib/python2.5/site-packages/scipy Python version 2.5.1 (r251:54863, Dec 21 2007, 09:21:07) [GCC 3.4.6 20060404 (Red Hat 3.4.6-3)] nose version 0.10.4 .................F........................................FF.F........................................................................................................................................................................................................................................................................................................... ====================================================================== FAIL: test_iv_cephes_vs_amos (test_basic.TestBessel)
Traceback (most recent call last): File "/data/home/nwagner/local/lib/python2.5/site-packages/scipy/special/tests/test_basic.py", line 1653, in test_iv_cephes_vs_amos self.check_cephes_vs_amos(iv, iv, rtol=1e-8, atol=1e-305) File "/data/home/nwagner/local/lib/python2.5/site-packages/scipy/special/tests/test_basic.py", line 1642, in check_cephes_vs_amos assert_tol_equal(c1, c2, err_msg=(v, z), rtol=rtol, atol=atol) File "/data/home/nwagner/local/lib/python2.5/site-packages/scipy/special/tests/test_basic.py", line 38, in assert_tol_equal verbose=verbose, header=header) File "/data/home/nwagner/local/lib/python2.5/site-packages/numpy/testing/utils.py", line 295, in assert_array_compare raise AssertionError(msg) AssertionError: Not equal to tolerance rtol=1e-08, atol=1e-305 (-120, -11) (mismatch 100.0%) x: array(1.3384173609003782e-110) y: array((1.3384173859242368e-110+0j)) ====================================================================== FAIL: test_yn_zeros (test_basic.TestBessel) ---------------------------------------------------------------------- Traceback (most recent call last): File "/data/home/nwagner/local/lib/python2.5/site-packages/scipy/special/tests/test_basic.py", line 1598, in test_yn_zeros 488.98055964441374646], rtol=1e-19) File "/data/home/nwagner/local/lib/python2.5/site-packages/scipy/special/tests/test_basic.py", line 38, in assert_tol_equal verbose=verbose, header=header) File "/data/home/nwagner/local/lib/python2.5/site-packages/numpy/testing/utils.py", line 295, in assert_array_compare raise AssertionError(msg) AssertionError: Not equal to tolerance rtol=1e-19, atol=0 (mismatch 100.0%) x: array([ 450.13573091, 463.05692375, 472.80651543, 481.27353185, 488.9805596 ]) y: array([ 450.13573092, 463.05692377, 472.80651546, 481.27353185, 488.98055964]) ====================================================================== FAIL: test_ynp_zeros (test_basic.TestBessel) ---------------------------------------------------------------------- Traceback (most recent call last): File "/data/home/nwagner/local/lib/python2.5/site-packages/scipy/special/tests/test_basic.py", line 1604, in test_ynp_zeros assert_tol_equal(yvp(443, ao), 0, atol=1e-15) File "/data/home/nwagner/local/lib/python2.5/site-packages/scipy/special/tests/test_basic.py", line 38, in assert_tol_equal verbose=verbose, header=header) File "/data/home/nwagner/local/lib/python2.5/site-packages/numpy/testing/utils.py", line 295, in assert_array_compare raise AssertionError(msg) AssertionError: Not equal to tolerance rtol=1e-07, atol=1e-15 (mismatch 100.0%) x: array([ 1.23909633e-10, -8.11850587e-16, 3.60822483e-16, 5.89805982e-16, 1.22644950e-15]) y: array(0) ====================================================================== FAIL: test_yv_cephes_vs_amos (test_basic.TestBessel) ---------------------------------------------------------------------- Traceback (most recent call last): File "/data/home/nwagner/local/lib/python2.5/site-packages/scipy/special/tests/test_basic.py", line 1650, in test_yv_cephes_vs_amos self.check_cephes_vs_amos(yv, yn, rtol=1e-11, atol=1e-305) File "/data/home/nwagner/local/lib/python2.5/site-packages/scipy/special/tests/test_basic.py", line 1640, in check_cephes_vs_amos assert c2.imag != 0, (v, z) AssertionError: (301, 1.0) ---------------------------------------------------------------------- Ran 361 tests in 0.782s FAILED (failures=4) <nose.result.TextTestResult run=361 errors=0 failures=4>
scipy.__version__ '0.8.0.dev5525'
Nils
Mon, 26 Jan 2009 11:27:30 +0100, Nils Wagner wrote:
====================================================================== FAIL: test_yv_cephes_vs_amos (test_basic.TestBessel) ---------------------------------------------------------------------- [clip] Traceback (most recent call last): File "/data/home/nwagner/local/lib/python2.5/site-packages/scipy/special/ tests/test_basic.py", line 1650, in test_yv_cephes_vs_amos self.check_cephes_vs_amos(yv, yn, rtol=1e-11, atol=1e-305) File "/data/home/nwagner/local/lib/python2.5/site-packages/scipy/special/ tests/test_basic.py", line 1640, in check_cephes_vs_amos assert c2.imag != 0, (v, z) AssertionError: (301, 1.0)
Appears to be a bug in cephes/yn on 64-bit systems:
import scipy import scipy.special as sc scipy.__version__ '0.7.0.dev5268' sc.yn(301,1) nan
whereas on a 32-bit system it gives the correct answer:
sc.yn(301,1) -inf
Pauli Virtanen wrote:
Mon, 26 Jan 2009 11:27:30 +0100, Nils Wagner wrote:
====================================================================== FAIL: test_yv_cephes_vs_amos (test_basic.TestBessel) ----------------------------------------------------------------------
[clip]
Traceback (most recent call last): File "/data/home/nwagner/local/lib/python2.5/site-packages/scipy/special/
tests/test_basic.py",
line 1650, in test_yv_cephes_vs_amos self.check_cephes_vs_amos(yv, yn, rtol=1e-11, atol=1e-305) File "/data/home/nwagner/local/lib/python2.5/site-packages/scipy/special/
tests/test_basic.py",
line 1640, in check_cephes_vs_amos assert c2.imag != 0, (v, z) AssertionError: (301, 1.0)
Appears to be a bug in cephes/yn on 64-bit systems:
import scipy import scipy.special as sc scipy.__version__
'0.7.0.dev5268'
sc.yn(301,1)
nan
whereas on a 32-bit system it gives the correct answer:
I wonder whether it would make sense to progressively use our own math function library, at least for some core routines - R core for example, does not use cephes or any external library (some functions are "translations" from specfun though). Or is there a simpler solution toward robust math routines (robustness against different architectures and out of domain input, in particular) ? David
Some of the real-valued Bessel function implementations from the Cephes library currently used in scipy.special have problems. (See #503, #851, #853, #854.) Fixing some of these (eg. #503) would require implementing robust computation algorithms from scratch. (The Specfun code is IMHO too obscure and badly commented to be relied on as an alternative.) However, the Boost library seems to have good implementations Bessel (and some other) special functions: http://svn.boost.org/svn/boost/trunk/boost/math/special_functions/detail/ http://www.boost.org/doc/libs/1_37_0/libs/math/doc/sf_and_dist/html/math_too... Also the license seems Scipy-compatible: http://www.boost.org/LICENSE_1_0.txt So, I'd like to bring these over to Scipy, to replace some of the Cephes routines. The only problem is that being in Boost, they are written in C++, and I guess we can't make Scipy to depend on it. I see two options: A) Bundle the relevant subset of Boost with Scipy. The problem here is that the special functions seem to pull in a sizable subset of the whole Boost library. Also, I don't know how well compilers handle the template-happy C++ in boost today on all platforms where Scipy must work on. B) Convert the Boost code from C++ to C. This is in fact quite trivial search-and-replace operation. One example here: http://github.com/pv/scipy/blob/ticket-503-special-iv-fix/scipy/special/ceph... I'd like to see (B) happen in scipy.special. Thoughts? -- Pauli Virtanen
The only problem is that being in Boost, they are written in C++, and I guess we can't make Scipy to depend on it.
The sparse module is already based on C+, so why not more ? Matthieu -- Information System Engineer, Ph.D. Website: http://matthieu-brucher.developpez.com/ Blogs: http://matt.eifelle.com and http://blog.developpez.com/?blog=92 LinkedIn: http://www.linkedin.com/in/matthieubrucher
Sun, 08 Feb 2009 18:25:32 +0100, Matthieu Brucher wrote:
The only problem is that being in Boost, they are written in C++, and I guess we can't make Scipy to depend on it.
The sparse module is already based on C+, so why not more ?
The problem is not C++ per se, but Boost: (i) How much of it we need to bundle with Scipy? (ii) Are there portability/build issues? But yes, using unmodified upstream code could be a relief from the maintenance POV. -- Pauli Virtanen
Sun, 08 Feb 2009 17:35:41 +0000, Pauli Virtanen wrote:
Sun, 08 Feb 2009 18:25:32 +0100, Matthieu Brucher wrote:
The only problem is that being in Boost, they are written in C++, and I guess we can't make Scipy to depend on it.
The sparse module is already based on C+, so why not more ?
The problem is not C++ per se, but Boost:
(i) How much of it we need to bundle with Scipy?
To answer myself: the bessel.hpp appears to pull in 2.7 Mb of boost headers. -- Pauli Virtanen
Pauli Virtanen wrote: Hi,
Sun, 08 Feb 2009 18:25:32 +0100, Matthieu Brucher wrote:
The only problem is that being in Boost, they are written in C++, and I guess we can't make Scipy to depend on it. The sparse module is already based on C+, so why not more ?
The problem is not C++ per se, but Boost:
(i) How much of it we need to bundle with Scipy? (ii) Are there portability/build issues?
Boost is a neverending souce of portability/build issues and every project I ever touched using boost had specific version requirements, i.e PuCUDA wanted either one of two speccfic release while quantlib wanted another set, but in between them there wasn't any boost that worked for both of them. Putting that code in-tree opens you up to all kinds of version mismatches and confusion if boost is installed system wide. I have had to fix or work around issues with recent boost on common platforms like OSX, much less seemingly "exotic" things like FreeBSD :), boost has its own build system (jam) which isn't exactly used commonly anywhere else and quite painful, i.e. boost always used the global Python headers for quantlib for example and you need either the latest release or some snapshot to work around that bug. Boost code requires beefy resources to compile and on and on an on. Please do not touch boost code, but if you must either translate C code or look at some alternative like mpmath, i.e. http://code.google.com/p/mpmath/ Fredrick is quite responsive about bugs and feature requests and we have talked to him about replacing some of the functionality provided by cephes in Sage via mpmath since they are arbitrary precision and pretty fast when optionally using gmp. But it also works in pure mode, i.e. all BSD licensed code.
But yes, using unmodified upstream code could be a relief from the maintenance POV.
Well, I am not so sure about that :) Cheers, Michael
Hi, On Sun, Feb 8, 2009 at 10:58 AM, Michael Abshoff <michael.abshoff@googlemail.com> wrote:
Pauli Virtanen wrote:
Hi,
Sun, 08 Feb 2009 18:25:32 +0100, Matthieu Brucher wrote:
The only problem is that being in Boost, they are written in C++, and I guess we can't make Scipy to depend on it. The sparse module is already based on C+, so why not more ?
The problem is not C++ per se, but Boost:
(i) How much of it we need to bundle with Scipy? (ii) Are there portability/build issues?
Boost is a neverending souce of portability/build issues and every project I ever touched using boost had specific version requirements, i.e PuCUDA wanted either one of two speccfic release while quantlib wanted another set, but in between them there wasn't any boost that worked for both of them. Putting that code in-tree opens you up to all kinds of version mismatches and confusion if boost is installed system wide.
I have had to fix or work around issues with recent boost on common platforms like OSX, much less seemingly "exotic" things like FreeBSD :), boost has its own build system (jam) which isn't exactly used commonly anywhere else and quite painful, i.e. boost always used the global Python headers for quantlib for example and you need either the latest release or some snapshot to work around that bug. Boost code requires beefy resources to compile and on and on an on. Please do not touch boost code, but if you must either translate C code or look at some alternative like mpmath, i.e.
http://code.google.com/p/mpmath/
Fredrick is quite responsive about bugs and feature requests and we have talked to him about replacing some of the functionality provided by cephes in Sage via mpmath since they are arbitrary precision and pretty fast when optionally using gmp. But it also works in pure mode, i.e. all BSD licensed code.
I completely agree with Michael here. Why not to use mpmath? It's bsd, it started as part of sympy and it was the GSoC project for sympy the last year. It's pretty competitive with gmp (e.g. for the pi digits calculations, it's even faster than Sage, unless Sage fixed that already), but one doesn't have to use gmp, if one doesn't want to. And I think both Fredrik and other mpmath and sympy developers would help to make mpmath working with scipy. Definitely I would. I think that's a better option, than to port some boost stuff and then you would have to maintain it. If you use mpmath, all of us win, imho. Ondrej
On Sun, Feb 8, 2009 at 18:12, Ondrej Certik <ondrej@certik.cz> wrote:
Hi,
On Sun, Feb 8, 2009 at 10:58 AM, Michael Abshoff <michael.abshoff@googlemail.com> wrote:
Pauli Virtanen wrote:
Hi,
Sun, 08 Feb 2009 18:25:32 +0100, Matthieu Brucher wrote:
The only problem is that being in Boost, they are written in C++, and I guess we can't make Scipy to depend on it. The sparse module is already based on C+, so why not more ?
The problem is not C++ per se, but Boost:
(i) How much of it we need to bundle with Scipy? (ii) Are there portability/build issues?
Boost is a neverending souce of portability/build issues and every project I ever touched using boost had specific version requirements, i.e PuCUDA wanted either one of two speccfic release while quantlib wanted another set, but in between them there wasn't any boost that worked for both of them. Putting that code in-tree opens you up to all kinds of version mismatches and confusion if boost is installed system wide.
I have had to fix or work around issues with recent boost on common platforms like OSX, much less seemingly "exotic" things like FreeBSD :), boost has its own build system (jam) which isn't exactly used commonly anywhere else and quite painful, i.e. boost always used the global Python headers for quantlib for example and you need either the latest release or some snapshot to work around that bug. Boost code requires beefy resources to compile and on and on an on. Please do not touch boost code, but if you must either translate C code or look at some alternative like mpmath, i.e.
http://code.google.com/p/mpmath/
Fredrick is quite responsive about bugs and feature requests and we have talked to him about replacing some of the functionality provided by cephes in Sage via mpmath since they are arbitrary precision and pretty fast when optionally using gmp. But it also works in pure mode, i.e. all BSD licensed code.
I completely agree with Michael here. Why not to use mpmath? It's bsd, it started as part of sympy and it was the GSoC project for sympy the last year. It's pretty competitive with gmp (e.g. for the pi digits calculations, it's even faster than Sage, unless Sage fixed that already), but one doesn't have to use gmp, if one doesn't want to.
And I think both Fredrik and other mpmath and sympy developers would help to make mpmath working with scipy. Definitely I would. I think that's a better option, than to port some boost stuff and then you would have to maintain it. If you use mpmath, all of us win, imho.
??? For implementing a C ufunc? -- Robert Kern "I have come to believe that the whole world is an enigma, a harmless enigma that is made terrible by our own mad attempt to interpret it as though it had an underlying truth." -- Umberto Eco
I completely agree with Michael here. Why not to use mpmath? It's bsd, it started as part of sympy and it was the GSoC project for sympy the last year. It's pretty competitive with gmp (e.g. for the pi digits calculations, it's even faster than Sage, unless Sage fixed that already), but one doesn't have to use gmp, if one doesn't want to.
And I think both Fredrik and other mpmath and sympy developers would help to make mpmath working with scipy. Definitely I would. I think that's a better option, than to port some boost stuff and then you would have to maintain it. If you use mpmath, all of us win, imho.
??? For implementing a C ufunc?
Using Cython? If it's as fast as anything else, why not. If it's not as fast, then that would be a reason not to use it. Ondrej
On Sun, Feb 8, 2009 at 21:38, Ondrej Certik <ondrej@certik.cz> wrote:
I completely agree with Michael here. Why not to use mpmath? It's bsd, it started as part of sympy and it was the GSoC project for sympy the last year. It's pretty competitive with gmp (e.g. for the pi digits calculations, it's even faster than Sage, unless Sage fixed that already), but one doesn't have to use gmp, if one doesn't want to.
And I think both Fredrik and other mpmath and sympy developers would help to make mpmath working with scipy. Definitely I would. I think that's a better option, than to port some boost stuff and then you would have to maintain it. If you use mpmath, all of us win, imho.
??? For implementing a C ufunc?
Using Cython? If it's as fast as anything else, why not. If it's not as fast, then that would be a reason not to use it.
Well, show me the code. :-) I wasn't aware that mpmath could be Cythonized. If it can, and the result is reasonably fast, that would be *really* useful. -- Robert Kern "I have come to believe that the whole world is an enigma, a harmless enigma that is made terrible by our own mad attempt to interpret it as though it had an underlying truth." -- Umberto Eco
Robert Kern wrote: Hi,
On Sun, Feb 8, 2009 at 21:38, Ondrej Certik <ondrej@certik.cz> wrote:
I completely agree with Michael here. Why not to use mpmath? It's bsd, it started as part of sympy and it was the GSoC project for sympy the last year. It's pretty competitive with gmp (e.g. for the pi digits calculations, it's even faster than Sage, unless Sage fixed that already), but one doesn't have to use gmp, if one doesn't want to.
And I think both Fredrik and other mpmath and sympy developers would help to make mpmath working with scipy. Definitely I would. I think that's a better option, than to port some boost stuff and then you would have to maintain it. If you use mpmath, all of us win, imho. ??? For implementing a C ufunc? Using Cython? If it's as fast as anything else, why not. If it's not as fast, then that would be a reason not to use it.
Well, show me the code. :-)
I wasn't aware that mpmath could be Cythonized. If it can, and the result is reasonably fast, that would be *really* useful.
I believe Ondrej was talking about using Cython to reduce call overhead, not to Cythonize mpmath which might or might not pay off. Given that Cepehes does not do arbitrary precision taking a look at mpmath before you decide to reinvent the wheel seems like a good idea. I don't really see the problem since mpmath works and in some cases is competitive with MPFR. And I don't meant that silly Pi to some 10^X computation which isn't particularly useful in the real world. It is quite hard to do arbitrary precision arithmetic and numerically stable special functions, so building on top of mpmath has its advantages. Cheers, Michael
??? For implementing a C ufunc? Using Cython? If it's as fast as anything else, why not. If it's not as fast, then that would be a reason not to use it.
Well, show me the code. :-)
I say that usually. :)
I wasn't aware that mpmath could be Cythonized. If it can, and the result is reasonably fast, that would be *really* useful.
I believe Ondrej was talking about using Cython to reduce call overhead, not to Cythonize mpmath which might or might not pay off.
Yes.
Given that Cepehes does not do arbitrary precision taking a look at mpmath before you decide to reinvent the wheel seems like a good idea. I don't really see the problem since mpmath works and in some cases is competitive with MPFR. And I don't meant that silly Pi to some 10^X computation which isn't particularly useful in the real world. It is quite hard to do arbitrary precision arithmetic and numerically stable special functions, so building on top of mpmath has its advantages.
Plus given that there are people who can help out with this. Ondrej
On Sun, Feb 8, 2009 at 22:30, Michael Abshoff <michael.abshoff@googlemail.com> wrote:
Robert Kern wrote:
Hi,
On Sun, Feb 8, 2009 at 21:38, Ondrej Certik <ondrej@certik.cz> wrote:
I completely agree with Michael here. Why not to use mpmath? It's bsd, it started as part of sympy and it was the GSoC project for sympy the last year. It's pretty competitive with gmp (e.g. for the pi digits calculations, it's even faster than Sage, unless Sage fixed that already), but one doesn't have to use gmp, if one doesn't want to.
And I think both Fredrik and other mpmath and sympy developers would help to make mpmath working with scipy. Definitely I would. I think that's a better option, than to port some boost stuff and then you would have to maintain it. If you use mpmath, all of us win, imho. ??? For implementing a C ufunc? Using Cython? If it's as fast as anything else, why not. If it's not as fast, then that would be a reason not to use it.
Well, show me the code. :-)
I wasn't aware that mpmath could be Cythonized. If it can, and the result is reasonably fast, that would be *really* useful.
I believe Ondrej was talking about using Cython to reduce call overhead, not to Cythonize mpmath which might or might not pay off.
Actually, that's all I meant by "Cythonize". You're right that it's not quite the right word.
Given that Cepehes does not do arbitrary precision taking a look at mpmath before you decide to reinvent the wheel seems like a good idea.
Actually, we weren't talking about multiprecision at all until you brought it up. The point of using the boost code wasn't the multiprecision aspect, but just that it was an alternative double-precision implementation that looked like it didn't have the bugs Cephes has. Multiprecision is an entirely different discussion, which goes some way towards explaining why I was confused why you and Ondrej think it's a good fit. -- Robert Kern "I have come to believe that the whole world is an enigma, a harmless enigma that is made terrible by our own mad attempt to interpret it as though it had an underlying truth." -- Umberto Eco
Robert Kern wrote: <SNIP>
I believe Ondrej was talking about using Cython to reduce call overhead, not to Cythonize mpmath which might or might not pay off.
Actually, that's all I meant by "Cythonize". You're right that it's not quite the right word.
Given that Cepehes does not do arbitrary precision taking a look at mpmath before you decide to reinvent the wheel seems like a good idea.
Actually, we weren't talking about multiprecision at all until you brought it up.
Fair enough. For the Sage people the fact that Cephes is limited in prevision and quite buggy is a major issue. That is why we are looking at mpmath.
The point of using the boost code wasn't the multiprecision aspect, but just that it was an alternative double-precision implementation that looked like it didn't have the bugs Cephes has.
By the way: In Sage we disabled all the special mtune and sse flags being set for the gfortran since it produced crashes for special functions all over the map on common architectures, i.e. Linux/P4 with gfortran 4.2.x as well as gfortran 4.3.x. On Solaris 10 I ended up using gfortran 4.3.2 for scipy and it works well, i.e. no crahes while gfortran 4.2.4 and earlier was a complete and buggy disaster. So someone might want to put a big warning in the release notes or the wiki that gfortran 4.2.x and Solaris do not play well together with numpy as well as scipy.
Multiprecision is an entirely different discussion, which goes some way towards explaining why I was confused why you and Ondrej think it's a good fit.
Those pesky non-engineering mathematicians and physicists :) Cheers, Michael
Michael Abshoff wrote:
By the way: In Sage we disabled all the special mtune and sse flags being set for the gfortran since it produced crashes for special functions all over the map on common architectures, i.e. Linux/P4 with gfortran 4.2.x as well as gfortran 4.3.x.
We disabled them as well in numpy trunk, so starting at 1.3, this should not be a problem. David
On Sun, Feb 8, 2009 at 10:23 AM, Pauli Virtanen <pav@iki.fi> wrote:
Some of the real-valued Bessel function implementations from the Cephes library currently used in scipy.special have problems. (See #503, #851, #853, #854.) Fixing some of these (eg. #503) would require implementing robust computation algorithms from scratch. (The Specfun code is IMHO too obscure and badly commented to be relied on as an alternative.)
However, the Boost library seems to have good implementations Bessel (and some other) special functions:
http://svn.boost.org/svn/boost/trunk/boost/math/special_functions/detail/
http://www.boost.org/doc/libs/1_37_0/libs/math/doc/sf_and_dist/html/math_too...
Also the license seems Scipy-compatible:
http://www.boost.org/LICENSE_1_0.txt
So, I'd like to bring these over to Scipy, to replace some of the Cephes routines.
The only problem is that being in Boost, they are written in C++, and I guess we can't make Scipy to depend on it.
I see two options:
A) Bundle the relevant subset of Boost with Scipy. The problem here is that the special functions seem to pull in a sizable subset of the whole Boost library.
That's a common problem with big c++ libraries.
Also, I don't know how well compilers handle the template-happy C++ in boost today on all platforms where Scipy must work on.
B) Convert the Boost code from C++ to C. This is in fact quite trivial search-and-replace operation. One example here:
http://github.com/pv/scipy/blob/ticket-503-special-iv-fix/scipy/special/ceph...
I'd like to see (B) happen in scipy.special. Thoughts?
I think it's a good idea. It would also be nice if we picked best of breed from several libraries to make up our own special functions collection. For instance, there are several log gamma functions. That's probably a big job though and we would need extensive tests for the functions before trying it. Do you know of any other project that has put together such a test suite? Does boost have versions for log1p? We need a better implementation in numpy itself. Chuck
Sun, 08 Feb 2009 11:39:55 -0700, Charles R Harris wrote: [clip]
I think it's a good idea. It would also be nice if we picked best of breed from several libraries to make up our own special functions collection. For instance, there are several log gamma functions.
I think we need to do this eventually, even if it means lots of work. At points the Cephes and Specfun codes seem like the author has not wanted to bother with the best possible algorithm, which leads to problems in corner cases.
That's probably a big job though and we would need extensive tests for the functions before trying it. Do you know of any other project that has put together such a test suite?
For Bessel functions we can easily test against the AMOS library, which appears to be reliable --- unless the order is close to negative integers in which case there can be cancellation errors of the order of 1e-6 in the reflection formulas. Boost itself has tests for its special functions, these are spot tests at precomputed points, from 1e-2..1e2 magnitudes in both parameters. [As an aside, I noticed Boost's I(v,x) overflows to infty somewhat earlier than necessary for very large orders, though.] GSL has similar point tests. (But it's GPLed.) Netlib/Specfun has a test suite: http://netlib.org/specfun/ ; in F77. Anyway, point tests across some magnitudes of parameters should be easy to generate. What takes more work is checking the behavior of the functions in transition regions where the method of computation changes, and asymptotic behavior (overflows, etc.) at large or small parameters and near singularities.
Does boost have versions for log1p? We need a better implementation in numpy itself.
It has. http://svn.boost.org/svn/boost/trunk/boost/math/special_functions/ log1p.hpp It's cluttered by C++ templates, but the algorithm looks like some serious effort has been put into it. -- Pauli Virtanen
On Mon, Feb 9, 2009 at 4:11 AM, Pauli Virtanen <pav@iki.fi> wrote:
Sun, 08 Feb 2009 11:39:55 -0700, Charles R Harris wrote: [clip]
I think it's a good idea. It would also be nice if we picked best of breed from several libraries to make up our own special functions collection. For instance, there are several log gamma functions.
I think we need to do this eventually, even if it means lots of work. At points the Cephes and Specfun codes seem like the author has not wanted to bother with the best possible algorithm, which leads to problems in corner cases.
Yes, I agree - I already asked about this a few weeks ago after some problems with other functions. I think cephes and specfun are not reliable - R does not use it, they have their own implementation of core math functions (sometimes inspired from cephes/specfun, but not that often).
Anyway, point tests across some magnitudes of parameters should be easy to generate. What takes more work is checking the behavior of the functions in transition regions where the method of computation changes, and asymptotic behavior (overflows, etc.) at large or small parameters and near singularities.
Yes, it would be a lot of work - I think we should focus on the tests before rewriting some functions. I would like to have a core scipy.special which is reliable: bessel, gamma/digamma/co, chebychev, basically most functions in R core would be a good start - and already quite heavy from a work POV.
It has.
http://svn.boost.org/svn/boost/trunk/boost/math/special_functions/ log1p.hpp
It's cluttered by C++ templates, but the algorithm looks like some serious effort has been put into it.
we could use their test-suite, maybe. cheers, David
On Feb 8, 2009, at 12:23 PM, Pauli Virtanen wrote:
However, the Boost library seems to have good implementations Bessel (and some other) special functions:
http://svn.boost.org/svn/boost/trunk/boost/math/special_functions/detail/ http://www.boost.org/doc/libs/1_37_0/libs/math/doc/sf_and_dist/html/math_too...
Also the license seems Scipy-compatible:
http://www.boost.org/LICENSE_1_0.txt
So, I'd like to bring these over to Scipy, to replace some of the Cephes routines.
The only problem is that being in Boost, they are written in C++, and I guess we can't make Scipy to depend on it.
What about the GSL implementation of the Bessel function? That's already in C and seems, in some sense, a more natural companion library to SciPy than boost. -gideon
On Sun, Feb 8, 2009 at 1:05 PM, Gideon Simpson <simpson@math.toronto.edu>wrote:
On Feb 8, 2009, at 12:23 PM, Pauli Virtanen wrote:
However, the Boost library seems to have good implementations Bessel (and some other) special functions:
http://svn.boost.org/svn/boost/trunk/boost/math/special_functions/detail/
http://www.boost.org/doc/libs/1_37_0/libs/math/doc/sf_and_dist/html/math_too...
Also the license seems Scipy-compatible:
http://www.boost.org/LICENSE_1_0.txt
So, I'd like to bring these over to Scipy, to replace some of the Cephes routines.
The only problem is that being in Boost, they are written in C++, and I guess we can't make Scipy to depend on it.
What about the GSL implementation of the Bessel function? That's already in C and seems, in some sense, a more natural companion library to SciPy than boost.
Wrong license. Chuck
On Feb 8, 2009, at 3:15 PM, Charles R Harris wrote:
What about the GSL implementation of the Bessel function? That's already in C and seems, in some sense, a more natural companion library to SciPy than boost.
Wrong license.
Chuck
I gather there is interest in reconciling the GPL and BSD licenses? -gideon
Gideon Simpson wrote:
On Feb 8, 2009, at 3:15 PM, Charles R Harris wrote:
What about the GSL implementation of the Bessel function? That's already in C and seems, in some sense, a more natural companion library to SciPy than boost.
Wrong license.
Chuck
I gather there is interest in reconciling the GPL and BSD licenses?
What do you mean by that? Without going into details and triggering an epic flamewar GPL and BSD define "freedom" differently and given 25 years of history between the FSF and the BSD camp and endless flamewars I don't see how this can be resolved in any way, shape or form. :)
-gideon
Cheers, Michael
_______________________________________________ Scipy-dev mailing list Scipy-dev@scipy.org http://projects.scipy.org/mailman/listinfo/scipy-dev
On Feb 8, 2009, at 3:42 PM, Michael Abshoff wrote:
What do you mean by that? Without going into details and triggering an epic flamewar GPL and BSD define "freedom" differently and given 25 years of history between the FSF and the BSD camp and endless flamewars I don't see how this can be resolved in any way, shape or form. :)
That's a real shame for an end-user like myself. I use a mixture of GSL and SciPy in my work, for different things, depending on the problem at hand. It would be good if there were robust interoperability there. -gideon
On Sun, Feb 08, 2009 at 03:58:45PM -0500, Gideon Simpson wrote:
On Feb 8, 2009, at 3:42 PM, Michael Abshoff wrote:
What do you mean by that? Without going into details and triggering an epic flamewar GPL and BSD define "freedom" differently and given 25 years of history between the FSF and the BSD camp and endless flamewars I don't see how this can be resolved in any way, shape or form. :)
That's a real shame for an end-user like myself. I use a mixture of GSL and SciPy in my work, for different things, depending on the problem at hand. It would be good if there were robust interoperability there.
That's because you don't have to worry about distributing software built upon these tools or making a profit. If you are trying to run a company, or simply if you are in a lab that wishes to sell some of the software it has developed, you end up caring about these things. Gaël
On Sun, Feb 8, 2009 at 14:58, Gideon Simpson <simpson@math.toronto.edu> wrote:
On Feb 8, 2009, at 3:42 PM, Michael Abshoff wrote:
What do you mean by that? Without going into details and triggering an epic flamewar GPL and BSD define "freedom" differently and given 25 years of history between the FSF and the BSD camp and endless flamewars I don't see how this can be resolved in any way, shape or form. :)
That's a real shame for an end-user like myself. I use a mixture of GSL and SciPy in my work, for different things, depending on the problem at hand. It would be good if there were robust interoperability there.
Interoperability is not the concern here. *You* can combine scipy and GSL all you like. That's not a problem. However, because we want to continue to use a BSD license for scipy, we don't include GPLed code like the GSL. -- Robert Kern "I have come to believe that the whole world is an enigma, a harmless enigma that is made terrible by our own mad attempt to interpret it as though it had an underlying truth." -- Umberto Eco
Gideon Simpson wrote:
On Feb 8, 2009, at 12:23 PM, Pauli Virtanen wrote:
However, the Boost library seems to have good implementations Bessel (and some other) special functions:
http://svn.boost.org/svn/boost/trunk/boost/math/special_functions/detail/
http://www.boost.org/doc/libs/1_37_0/libs/math/doc/sf_and_dist/html/math_too...
Also the license seems Scipy-compatible:
http://www.boost.org/LICENSE_1_0.txt
So, I'd like to bring these over to Scipy, to replace some of the Cephes routines.
The only problem is that being in Boost, they are written in C++, and I guess we can't make Scipy to depend on it.
Hi,
What about the GSL implementation of the Bessel function? That's already in C and seems, in some sense, a more natural companion library to SciPy than boost.
GSL is GPL licensed - not surprisingly since it is the GNU scientifc library :)
-gideon
Cheers, Michael
------------------------------------------------------------------------
_______________________________________________ Scipy-dev mailing list Scipy-dev@scipy.org http://projects.scipy.org/mailman/listinfo/scipy-dev
On Mon, Feb 9, 2009 at 2:23 AM, Pauli Virtanen <pav@iki.fi> wrote:
Some of the real-valued Bessel function implementations from the Cephes library currently used in scipy.special have problems. (See #503, #851, #853, #854.) Fixing some of these (eg. #503) would require implementing robust computation algorithms from scratch. (The Specfun code is IMHO too obscure and badly commented to be relied on as an alternative.)
However, the Boost library seems to have good implementations Bessel (and some other) special functions:
http://svn.boost.org/svn/boost/trunk/boost/math/special_functions/detail/ http://www.boost.org/doc/libs/1_37_0/libs/math/doc/sf_and_dist/html/math_too...
Also the license seems Scipy-compatible:
http://www.boost.org/LICENSE_1_0.txt
So, I'd like to bring these over to Scipy, to replace some of the Cephes routines.
The only problem is that being in Boost, they are written in C++, and I guess we can't make Scipy to depend on it.
I see two options:
A) Bundle the relevant subset of Boost with Scipy. The problem here is that the special functions seem to pull in a sizable subset of the whole Boost library.
Also, I don't know how well compilers handle the template-happy C++ in boost today on all platforms where Scipy must work on.
I am -1 on boost. It is a nightmare to support on many platforms, and it is unreadable for people who are not C++ hackers.
B) Convert the Boost code from C++ to C. This is in fact quite trivial search-and-replace operation. One example here:
http://github.com/pv/scipy/blob/ticket-503-special-iv-fix/scipy/special/ceph...
I'd like to see (B) happen in scipy.special. Thoughts?
This is much better - I really don't see the point of using C++ for math functions. I am ok with this. David
David Cournapeau wrote: Hi David, <SNIP>
This is much better - I really don't see the point of using C++ for math functions. I am ok with this.
David
Out of curiosity: I checked the boost website and it states for the math lib: "All the implementations are fully generic and support the use of arbitrary "real-number" types, although they are optimised for use with types with known-about significand (or mantissa) sizes: typically float, double or long double." Since I assume some people around here are interested in arbitrary precisions and after looking some more at the documentation it seems that that library only supports this via using an NTL type which in turn uses GMP. NTL itself is GPLed, GMP is LGPL, so either one does not fit the licensing requirements of Scipy. David mentioned to write a library from scratch and I also assume that you want arbitrary precisions. Given that GMP is LGPL, the arbitrary precisions code in OpenSSL is covered by a BSD advertising clause (which might or might not be a deal breaker around here) what do you suggest to do about arbitrary precisions? I am not aware of any BSD 2 or 3 clause license library besides mpmath :) Cheers, Michael
_______________________________________________ Scipy-dev mailing list Scipy-dev@scipy.org http://projects.scipy.org/mailman/listinfo/scipy-dev
Hi Michael, Michael Abshoff wrote:
Since I assume some people around here are interested in arbitrary precisions and after looking some more at the documentation it seems that that library only supports this via using an NTL type which in turn uses GMP. NTL itself is GPLed, GMP is LGPL, so either one does not fit the licensing requirements of Scipy.
David mentioned to write a library from scratch and I also assume that you want arbitrary precisions.
Actually, I did not think about arbitrary precision at all. Not that it would not be nice, but I simply know nothing about the topic :) I am far from being entirely convinced by my own suggestion of writing from scratch: it is a huge amount of work, and the possibilities to get it wrong are numerous. But I am wondering whether we have a choice: the only existing code which is licence compatible with scipy is the one we use now (cephes, specfun, etc...), and the code is not great (cephes is not even ANSI C, for example), except for toms maybe. My suggestion is based on how R is doing things; I tend to consider R as a reference, or at least a pretty good baseline when precision and accuracy are concerned. R does not use Cephes nor specfun. We of course can't use R code, but we could at least take inspiration of their references (e.g. citations for implementation), and use R for testing. I also think one of the problem of scipy.special is its size - there are so many functions, with little to no testing. So maybe we could make up a list of a small subset of functions which are 'essential', and focus on them (we would of course keep the current code). I don't know whether such a small subset exists. For information, nmath (the core maths routines in R) is ~7500 LOC according to sloccount, and they have ~ 100 functions, maybe that would be a good subset ? I also have no idea how to test those functions: when someones says their function is precise up to 1e-6, is it against some theoretical values which are computable (like gamma(0.5) = euler constant), from theoretical consideration on the implementation ? Pauli talked about point tests, but it sounds hard to get the right points for testing ? cheers, David
Sun, 08 Feb 2009 20:25:37 -0800, Michael Abshoff wrote:
David Cournapeau wrote:
This is much better - I really don't see the point of using C++ for math functions. I am ok with this.
Out of curiosity: I checked the boost website and it states for the math lib:
"All the implementations are fully generic and support the use of arbitrary "real-number" types, although they are optimised for use with types with known-about significand (or mantissa) sizes: typically float, double or long double."
Since I assume some people around here are interested in arbitrary precisions and after looking some more at the documentation it seems that that library only supports this via using an NTL type which in turn uses GMP. NTL itself is GPLed, GMP is LGPL, so either one does not fit the licensing requirements of Scipy.
Yeh, arbitrary precision could be nice in principle. But as I see it, at the moment it's out-of-scope for Scipy. Right now, we only need good implementations in double precision. These we can get for some functions for example by adapting the Boost code (and re-testing it) -- this is much less work than rewriting everything from scratch.
David mentioned to write a library from scratch and I also assume that you want arbitrary precisions. Given that GMP is LGPL, the arbitrary precisions code in OpenSSL is covered by a BSD advertising clause (which might or might not be a deal breaker around here) what do you suggest to do about arbitrary precisions? I am not aware of any BSD 2 or 3 clause license library besides mpmath :)
For the present, I'd say that we should leave the arbitrary-precision functions implemented in mpmath. -- Pauli Virtanen
On Mon, Feb 9, 2009 at 12:35 AM, Pauli Virtanen <pav@iki.fi> wrote:
Sun, 08 Feb 2009 20:25:37 -0800, Michael Abshoff wrote:
David Cournapeau wrote:
This is much better - I really don't see the point of using C++ for math functions. I am ok with this.
Out of curiosity: I checked the boost website and it states for the math lib:
"All the implementations are fully generic and support the use of arbitrary "real-number" types, although they are optimised for use with types with known-about significand (or mantissa) sizes: typically float, double or long double."
Since I assume some people around here are interested in arbitrary precisions and after looking some more at the documentation it seems that that library only supports this via using an NTL type which in turn uses GMP. NTL itself is GPLed, GMP is LGPL, so either one does not fit the licensing requirements of Scipy.
Yeh, arbitrary precision could be nice in principle.
But as I see it, at the moment it's out-of-scope for Scipy. Right now, we only need good implementations in double precision. These we can get for some functions for example by adapting the Boost code (and re-testing it) -- this is much less work than rewriting everything from scratch.
David mentioned to write a library from scratch and I also assume that you want arbitrary precisions. Given that GMP is LGPL, the arbitrary precisions code in OpenSSL is covered by a BSD advertising clause (which might or might not be a deal breaker around here) what do you suggest to do about arbitrary precisions? I am not aware of any BSD 2 or 3 clause license library besides mpmath :)
For the present, I'd say that we should leave the arbitrary-precision functions implemented in mpmath.
Right. For double precision I think mpmath is not so fast. Fredrik, is it difficult to make mpmath fast even for double precision? Last time I asked: http://groups.google.com/group/mpmath/browse_thread/thread/bca53c3382945c34/ you replied: " SciPy already provides a truckload of machine precision special functions, with excellent (fast and robust) implementations. It'd be hard to top that. " But apparently, maybe mpmath can be useful. Ondrej
Mon, 09 Feb 2009 01:20:11 -0800, Ondrej Certik wrote: [clip]
Right. For double precision I think mpmath is not so fast. Fredrik, is it difficult to make mpmath fast even for double precision?
Last time I asked:
http://groups.google.com/group/mpmath/browse_thread/thread/ bca53c3382945c34/
you replied:
" SciPy already provides a truckload of machine precision special functions, with excellent (fast and robust) implementations. It'd be hard to top that. "
But apparently, maybe mpmath can be useful.
I'd say that mpmath faces the same robustness and testing issues as Scipy with regard to special functions. (In addition, since it's written in Python, I'd assume it also faces additional performance issues.) Also, algorithms that work well in arbitrary precision might not work for limited precision, due to loss of precision or under/overflows in intermediate steps. Looking at the Bessel function implementations in mpmath/functions.py, I'd say that at least besselj and besseli would face overflow issues for large arguments if they were working in double precision. This kind of issues are actually the most difficult to get right. To clarify: Definitely I think that mpmath is great work, and I'm happy to see people working on it, including improvements to its special function library. But at the present, I think the path of least resistance for Scipy is to continue using, testing, and improving existing implementations of special function codes, written in C or F77, directly for limited precision. -- Pauli Virtanen
On Mon, Feb 9, 2009 at 3:26 AM, Pauli Virtanen <pav@iki.fi> wrote:
Mon, 09 Feb 2009 01:20:11 -0800, Ondrej Certik wrote: [clip]
Right. For double precision I think mpmath is not so fast. Fredrik, is it difficult to make mpmath fast even for double precision?
Last time I asked:
http://groups.google.com/group/mpmath/browse_thread/thread/ bca53c3382945c34/<http://groups.google.com/group/mpmath/browse_thread/thread/%0Abca53c3382945c...>
you replied:
" SciPy already provides a truckload of machine precision special functions, with excellent (fast and robust) implementations. It'd be hard to top that. "
But apparently, maybe mpmath can be useful.
I'd say that mpmath faces the same robustness and testing issues as Scipy with regard to special functions. (In addition, since it's written in Python, I'd assume it also faces additional performance issues.)
Also, algorithms that work well in arbitrary precision might not work for limited precision, due to loss of precision or under/overflows in intermediate steps. Looking at the Bessel function implementations in mpmath/functions.py, I'd say that at least besselj and besseli would face overflow issues for large arguments if they were working in double precision. This kind of issues are actually the most difficult to get right.
To clarify: Definitely I think that mpmath is great work, and I'm happy to see people working on it, including improvements to its special function library.
But at the present, I think the path of least resistance for Scipy is to continue using, testing, and improving existing implementations of special function codes, written in C or F77, directly for limited precision.
Here are some references for both testing and implementations: http://math.nist.gov/mcsd/Reports/2001/nesf/ . I found Cody's book(s) a good reference back in the day. I wonder if anything happened with this proposal: http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.44.7298 ? Note the references to Cody which are to works with both implementations and tests. I think Cody was responsible for specfun, have there been problems with that package? Chuck
Mon, 09 Feb 2009 09:28:42 -0700, Charles R Harris wrote: [clip]
Here are some references for both testing and implementations: http://math.nist.gov/mcsd/Reports/2001/nesf/ . I found Cody's book(s) a good reference back in the day.
I wonder if anything happened with this proposal: http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.44.7298 ?
Thanks for the references.
Note the references to Cody which are to works with both implementations and tests. I think Cody was responsible for specfun, have there been problems with that package?
Confusingly enough, the 'specfun.f' in Scipy is AFAICS not Cody's netlib.org/specfun. Cody's code does have tests. -- Pauli Virtanen
On Mon, Feb 9, 2009 at 10:20 AM, Ondrej Certik <ondrej@certik.cz> wrote:
Right. For double precision I think mpmath is not so fast. Fredrik, is it difficult to make mpmath fast even for double precision?
Fast and accurate machine floating-point implementations require entirely different algorithms. One also needs to piece together multiple algorithms to handle cases where, in arbitrary-precision arithmetic, increasing some parameter is sufficient. I can see a use for a pure-Python double precision special functions library. I don't think mpmath should be it, at least not for now; getting the arbitrary-precision algorithms right is enough work. Should someone be interested in writing such a library, they could use mpmath to test accuracy against, and even to generate Chebyshev approximations and the like. Fredrik
Pauli Virtanen wrote:
Sun, 08 Feb 2009 20:25:37 -0800, Michael Abshoff wrote:
David Cournapeau wrote:
This is much better - I really don't see the point of using C++ for math functions. I am ok with this.
Out of curiosity: I checked the boost website and it states for the math lib:
"All the implementations are fully generic and support the use of arbitrary "real-number" types, although they are optimised for use with types with known-about significand (or mantissa) sizes: typically float, double or long double."
Since I assume some people around here are interested in arbitrary precisions and after looking some more at the documentation it seems that that library only supports this via using an NTL type which in turn uses GMP. NTL itself is GPLed, GMP is LGPL, so either one does not fit the licensing requirements of Scipy.
Yeh, arbitrary precision could be nice in principle.
But as I see it, at the moment it's out-of-scope for Scipy. Right now, we only need good implementations in double precision. These we can get for some functions for example by adapting the Boost code (and re-testing it) -- this is much less work than rewriting everything from scratch.
Are you familiar with boost testing ? I wonder whether it would be possible to automatically convert it to something usable for scipy (the .ipp files which contain the data should be relatively easy to convert to python, they all look the same with almost no code at all). cheers, David
On Mon, Feb 9, 2009 at 7:35 PM, David Cournapeau <david@ar.media.kyoto-u.ac.jp> wrote:
Pauli Virtanen wrote:
Sun, 08 Feb 2009 20:25:37 -0800, Michael Abshoff wrote:
David Cournapeau wrote:
This is much better - I really don't see the point of using C++ for math functions. I am ok with this.
Out of curiosity: I checked the boost website and it states for the math lib:
"All the implementations are fully generic and support the use of arbitrary "real-number" types, although they are optimised for use with types with known-about significand (or mantissa) sizes: typically float, double or long double."
Since I assume some people around here are interested in arbitrary precisions and after looking some more at the documentation it seems that that library only supports this via using an NTL type which in turn uses GMP. NTL itself is GPLed, GMP is LGPL, so either one does not fit the licensing requirements of Scipy.
Yeh, arbitrary precision could be nice in principle.
But as I see it, at the moment it's out-of-scope for Scipy. Right now, we only need good implementations in double precision. These we can get for some functions for example by adapting the Boost code (and re-testing it) -- this is much less work than rewriting everything from scratch.
Are you familiar with boost testing ? I wonder whether it would be possible to automatically convert it to something usable for scipy (the .ipp files which contain the data should be relatively easy to convert to python, they all look the same with almost no code at all).
I started to work on this, actually. Getting the test data is easy, but there is still some manual work to know which function is used for which data (it does not look like the files are consistent enough so that the mapping data file -> function tested can be done automatically). Shall I integrate this into scipy ? cheers, David
On Tue, Feb 10, 2009 at 5:39 AM, David Cournapeau <cournape@gmail.com>wrote:
On Mon, Feb 9, 2009 at 7:35 PM, David Cournapeau <david@ar.media.kyoto-u.ac.jp> wrote:
Pauli Virtanen wrote:
Sun, 08 Feb 2009 20:25:37 -0800, Michael Abshoff wrote:
David Cournapeau wrote:
This is much better - I really don't see the point of using C++ for math functions. I am ok with this.
Out of curiosity: I checked the boost website and it states for the math lib:
"All the implementations are fully generic and support the use of arbitrary "real-number" types, although they are optimised for use with types with known-about significand (or mantissa) sizes: typically float, double or long double."
Since I assume some people around here are interested in arbitrary precisions and after looking some more at the documentation it seems that that library only supports this via using an NTL type which in turn uses GMP. NTL itself is GPLed, GMP is LGPL, so either one does not fit the licensing requirements of Scipy.
Yeh, arbitrary precision could be nice in principle.
But as I see it, at the moment it's out-of-scope for Scipy. Right now, we only need good implementations in double precision. These we can get for some functions for example by adapting the Boost code (and re-testing it) -- this is much less work than rewriting everything from scratch.
Are you familiar with boost testing ? I wonder whether it would be possible to automatically convert it to something usable for scipy (the .ipp files which contain the data should be relatively easy to convert to python, they all look the same with almost no code at all).
I started to work on this, actually. Getting the test data is easy, but there is still some manual work to know which function is used for which data (it does not look like the files are consistent enough so that the mapping data file -> function tested can be done automatically).
Shall I integrate this into scipy ?
I think you should also take a look at Cody's netlib.org/specfun. It has the Bessel functions along with the tests and Cody does careful work with attention to detail and error. As Pauli says, the netlib specfunc isn't the same as the specfunc in Scipy. The notes in specfun are also worth a read. Chuck
Tue, 10 Feb 2009 21:39:09 +0900, David Cournapeau wrote: [clip: tests for special functions]
I started to work on this, actually. Getting the test data is easy, but there is still some manual work to know which function is used for which data (it does not look like the files are consistent enough so that the mapping data file -> function tested can be done automatically).
Shall I integrate this into scipy ?
Please do. It can't hurt to have better tests. And we can have them committed even before we start to fix any bugs they expose. Btw, it would be nice if buildbot.scipy.org also handled Scipy in addition to Numpy. (Buildbot unfortunately has no real support for multiple projects, so this would require running a second buildmaster daemon on a separate port.) Meanwhile, there's a dump of local buildbot results here: http://www.iki.fi/pav/tmp/bb/scipy/waterfall/ but this machine can't really be used as a build master as the machine it isn't online 24/7. -- Pauli Virtanen
On Wed, Feb 11, 2009 at 1:15 AM, Pauli Virtanen <pav@iki.fi> wrote:
Tue, 10 Feb 2009 21:39:09 +0900, David Cournapeau wrote: [clip: tests for special functions]
I started to work on this, actually. Getting the test data is easy, but there is still some manual work to know which function is used for which data (it does not look like the files are consistent enough so that the mapping data file -> function tested can be done automatically).
Shall I integrate this into scipy ?
Please do. It can't hurt to have better tests. And we can have them committed even before we start to fix any bugs they expose.
I started a branch, special_refactor. I added all the converted Boost data set (the .ipp files to .csv), plus the small python script I used to generate them. I started implementing the corresponding tests - but this takes some time, because of all this template stuff which is awkward to follow. The only thing to do is to find which function is called for which test with which parameter - someone more familiar with boost could to this much faster, I guess. David
Wed, 11 Feb 2009 03:31:30 +0900, David Cournapeau wrote: [clip]
I started a branch, special_refactor. I added all the converted Boost data set (the .ipp files to .csv), plus the small python script I used to generate them. I started implementing the corresponding tests - but this takes some time, because of all this template stuff which is awkward to follow. The only thing to do is to find which function is called for which test with which parameter - someone more familiar with boost could to this much faster, I guess.
I added a couple of more functions to the tests: They correctly point out that in 0.7.0: + The problems in Cephes's Iv (large argument), Yv (large order) and Kn (large order) + Numpy's complex-valued `arcsinh` and `arctanh` can have large relative errors (~1e-5) for small arguments (< eps)! Loss of precision in the naive implementation, I'll bet. but they fail to spot the other known issues. But on the positive side, the `arcsinh` issue is the only new one that came up. One problem with these tests is that the data files are *huge*, they currently total ~ 7 Mb. Even compressed, or saved as .npy files, these would add ~ 2 Mb to the Scipy source tarball. So I'm not sure what to do with this... -- Pauli Virtanen
Pauli Virtanen wrote:
Wed, 11 Feb 2009 03:31:30 +0900, David Cournapeau wrote:
[clip]
I started a branch, special_refactor. I added all the converted Boost data set (the .ipp files to .csv), plus the small python script I used to generate them. I started implementing the corresponding tests - but this takes some time, because of all this template stuff which is awkward to follow. The only thing to do is to find which function is called for which test with which parameter - someone more familiar with boost could to this much faster, I guess.
I added a couple of more functions to the tests:
They correctly point out that in 0.7.0:
+ The problems in Cephes's Iv (large argument), Yv (large order) and Kn (large order)
+ Numpy's complex-valued `arcsinh` and `arctanh` can have large relative errors (~1e-5) for small arguments (< eps)!
Loss of precision in the naive implementation, I'll bet.
but they fail to spot the other known issues. But on the positive side, the `arcsinh` issue is the only new one that came up.
One problem with these tests is that the data files are *huge*, they currently total ~ 7 Mb. Even compressed, or saved as .npy files, these would add ~ 2 Mb to the Scipy source tarball. So I'm not sure what to do with this...
That's the reason why I started a branch - I did not know how it would end up. I don't see an obvious answer to the problem: those are tests for ~ 100 functions, so this means 20kb of compressed data/function on average. Each test is two data points at least (x and f(x)), this means around ~ 500 test points/function. That does not sound that big anymore. Maybe we could have an option to split the dataset to make them separate from the main tarball ? I kept the data in .csv because I thought it would be nice to test for double and float at least, and the gain using binary would not be that huge anymore (it is also easier to use for tests outside the python machinery), cheers, David
On Thu, Feb 12, 2009 at 6:31 PM, David Cournapeau < david@ar.media.kyoto-u.ac.jp> wrote:
Pauli Virtanen wrote:
Wed, 11 Feb 2009 03:31:30 +0900, David Cournapeau wrote:
[clip]
I started a branch, special_refactor. I added all the converted Boost data set (the .ipp files to .csv), plus the small python script I used to generate them. I started implementing the corresponding tests - but this takes some time, because of all this template stuff which is awkward to follow. The only thing to do is to find which function is called for which test with which parameter - someone more familiar with boost could to this much faster, I guess.
I added a couple of more functions to the tests:
They correctly point out that in 0.7.0:
+ The problems in Cephes's Iv (large argument), Yv (large order) and Kn (large order)
+ Numpy's complex-valued `arcsinh` and `arctanh` can have large relative errors (~1e-5) for small arguments (< eps)!
Loss of precision in the naive implementation, I'll bet.
but they fail to spot the other known issues. But on the positive side, the `arcsinh` issue is the only new one that came up.
One problem with these tests is that the data files are *huge*, they currently total ~ 7 Mb. Even compressed, or saved as .npy files, these would add ~ 2 Mb to the Scipy source tarball. So I'm not sure what to do with this...
That's the reason why I started a branch - I did not know how it would end up. I don't see an obvious answer to the problem: those are tests for ~ 100 functions, so this means 20kb of compressed data/function on average. Each test is two data points at least (x and f(x)), this means around ~ 500 test points/function. That does not sound that big anymore. Maybe we could have an option to split the dataset to make them separate from the main tarball ?
I kept the data in .csv because I thought it would be nice to test for double and float at least, and the gain using binary would not be that huge anymore (it is also easier to use for tests outside the python machinery),
Maybe it would be best to split out the tests into a separate project and not distribute it with scipy. It could be turned into a generic test suite based on python that could be used to test any implementation of a specific function. Chuck
On Sun, 25 Jan 2009 20:03:26 +0100 "Nils Wagner" <nwagner@iam.uni-stuttgart.de> wrote:
Hi all,
I do get three failures wrt. to special functions and python2.6
====================================================================== FAIL: test_yn_zeros (test_basic.TestBessel) ---------------------------------------------------------------------- Traceback (most recent call last): File "/home/nwagner/local/lib64/python2.6/site-packages/scipy/special/tests/test_basic.py", line 1598, in test_yn_zeros 488.98055964441374646], rtol=1e-19) File "/home/nwagner/local/lib64/python2.6/site-packages/scipy/special/tests/test_basic.py", line 38, in assert_tol_equal verbose=verbose, header=header) File "/home/nwagner/local/lib64/python2.6/site-packages/numpy/testing/utils.py", line 295, in assert_array_compare raise AssertionError(msg) AssertionError: Not equal to tolerance rtol=1e-19, atol=0
(mismatch 100.0%) x: array([ 450.136, 463.057, 472.807, 481.274, 488.981]) y: array([ 450.136, 463.057, 472.807, 481.274, 488.981])
====================================================================== FAIL: test_ynp_zeros (test_basic.TestBessel) ---------------------------------------------------------------------- Traceback (most recent call last): File "/home/nwagner/local/lib64/python2.6/site-packages/scipy/special/tests/test_basic.py", line 1604, in test_ynp_zeros assert_tol_equal(yvp(443, ao), 0, atol=1e-15) File "/home/nwagner/local/lib64/python2.6/site-packages/scipy/special/tests/test_basic.py", line 38, in assert_tol_equal verbose=verbose, header=header) File "/home/nwagner/local/lib64/python2.6/site-packages/numpy/testing/utils.py", line 295, in assert_array_compare raise AssertionError(msg) AssertionError: Not equal to tolerance rtol=1e-07, atol=1e-15
(mismatch 100.0%) x: array([ 1.239e-10, -8.119e-16, 3.608e-16, 5.898e-16, 1.226e-15]) y: array(0)
====================================================================== FAIL: Negative-order Bessels ---------------------------------------------------------------------- Traceback (most recent call last): File "/home/nwagner/local/lib64/python2.6/site-packages/nose-0.10.4-py2.6.egg/nose/case.py", line 182, in runTest self.test(*self.arg) File "/home/nwagner/local/lib64/python2.6/site-packages/scipy/special/tests/test_basic.py", line 1665, in test_ticket_853 assert_tol_equal(iv(-0.5, 1 ), 1.231200214592967) File "/home/nwagner/local/lib64/python2.6/site-packages/scipy/special/tests/test_basic.py", line 38, in assert_tol_equal verbose=verbose, header=header) File "/home/nwagner/local/lib64/python2.6/site-packages/numpy/testing/utils.py", line 295, in assert_array_compare raise AssertionError(msg) AssertionError: Not equal to tolerance rtol=1e-07, atol=0
(mismatch 100.0%) x: array(0.0) y: array(1.231200214592967)
Nils _______________________________________________ Scipy-dev mailing list Scipy-dev@scipy.org http://projects.scipy.org/mailman/listinfo/scipy-dev
Another failure is present in
scipy.__version__ '0.8.0.dev5525'
====================================================================== FAIL: test_iv_cephes_vs_amos (test_basic.TestBessel) ---------------------------------------------------------------------- Traceback (most recent call last): File "/data/home/nwagner/local/lib/python2.5/site-packages/scipy/special/tests/test_basic.py", line 1653, in test_iv_cephes_vs_amos self.check_cephes_vs_amos(iv, iv, rtol=1e-8, atol=1e-305) File "/data/home/nwagner/local/lib/python2.5/site-packages/scipy/special/tests/test_basic.py", line 1642, in check_cephes_vs_amos assert_tol_equal(c1, c2, err_msg=(v, z), rtol=rtol, atol=atol) File "/data/home/nwagner/local/lib/python2.5/site-packages/scipy/special/tests/test_basic.py", line 38, in assert_tol_equal verbose=verbose, header=header) File "/data/home/nwagner/local/lib/python2.5/site-packages/numpy/testing/utils.py", line 295, in assert_array_compare raise AssertionError(msg) AssertionError: Not equal to tolerance rtol=1e-08, atol=1e-305 (-120, -11) (mismatch 100.0%) x: array(1.3384173609003782e-110) y: array((1.3384173859242368e-110+0j))
Mon, 26 Jan 2009 09:03:19 +0100, Nils Wagner wrote: [clip]
====================================================================== FAIL: test_iv_cephes_vs_amos (test_basic.TestBessel) ---------------------------------------------------------------------- [clip] Not equal to tolerance rtol=1e-08, atol=1e-305 (-120, -11) (mismatch 100.0%) x: array(1.3384173609003782e-110) y: array((1.3384173859242368e-110+0j))
The way cephes/iv computes the value of Iv is not too accurate for large orders or large x. I have some code to compute the result in this range from proper asymptotic expansions, but also it isn't committed yet. -- Pauli Virtanen
Sun, 25 Jan 2009 20:03:26 +0100, Nils Wagner wrote:
I do get three failures wrt. to special functions and python2.6
Too strict tolerances in the tests for the most part, I think. Will fix. I presume this is a 64-bit Linux platform. I'm in progress of addressing several bugs and improving the tests for the Bessel-related special functions in Scipy trunk. If someone wants to review the changesets I've committed recently, please go ahead.
====================================================================== FAIL: test_yn_zeros (test_basic.TestBessel) ---------------------------------------------------------------------- [clip] Not equal to tolerance rtol=1e-19, atol=0
Too strict tolerance, I think. Some errors ~eps probably arise, depending on the compiler.
====================================================================== FAIL: test_ynp_zeros (test_basic.TestBessel) ---------------------------------------------------------------------- [clip] AssertionError: Not equal to tolerance rtol=1e-07, atol=1e-15
(mismatch 100.0%) x: array([ 1.239e-10, -8.119e-16, 3.608e-16, 5.898e-16, 1.226e-15]) y: array(0)
Ditto. Except maybe for the first item. The question now is whether the problem is in the `yvp` routine or in `ynp_zeros`.
====================================================================== FAIL: Negative-order Bessels ---------------------------------------------------------------------- [clip] assert_tol_equal(iv(-0.5, 1 ), 1.231200214592967) [clip]
Negative half-integers are a known failure of cephes/iv (also Scipy 0.6.0), since hyperg has a pole at the corresponding points. But I think this test is commented out in the current SVN trunk, though. I have a fix for this, but it's not committed yet. -- Pauli Virtanen
participants (12)
-
Charles R Harris -
David Cournapeau -
David Cournapeau -
Fredrik Johansson -
Gael Varoquaux -
Gideon Simpson -
Matthieu Brucher -
Michael Abshoff -
Nils Wagner -
Ondrej Certik -
Pauli Virtanen -
Robert Kern