Suggestions for GSoC Projects
![](https://secure.gravatar.com/avatar/bae0c4adb1469b38cfe4b1844180c446.jpg?s=120&d=mm&r=g)
With GSoC 2014 being round the corner, I hereby put up few projects for discussion that I would love to pursue as a student. Guidance, suggestions are cordially welcome:- 1. If I am not mistaken, contour integration is not supported by SciPy; in fact even line integrals of real functions is yet to be implemented in SciPy, which is surprising. Though we at present have SymPy for line Integrals, I doubt if there is any open-source python package supporting the calculation of Contour Integrals. With integrate module of SciPy already having been properly developed for definite integration, implementation of line as well as contour integrals, I presume; would not require work from scratch and shall be a challenging but fruitful project. 2. I really have no idea if the purpose of NumPy or SciPy would encompass this but we are yet to have indefinite integration. An implementation of that, though highly challenging, may open doors for innumerable other functions like the ones to calculate the Laplace transform, Hankel transform and many more. 3. As stated earlier, we have spherical harmonic functions (with much scope for dev) we are yet to have elliptical and cylindrical harmonic function, which may be developed. 4. Lastly, we are yet to have Inverse Laplace transforms which as Ralf has rightly pointed out it may be too challenging to implement. 5. Further reading the road-map given by Mr.Ralf, I would like to develop the Bluestein's FFT algorithm. Thanks for reading along till the end. I shall append to this mail as when I am struck with ideas. Please do give your valuable guidance
![](https://secure.gravatar.com/avatar/af6c39d6943bd4b0e1fde23161e7bb8c.jpg?s=120&d=mm&r=g)
On Fri, 31 Jan 2014 04:31:01 +0530, jennifer stone wrote:
As stated before, I am personally interested in seeing the spherical harmonics in SciPy improve.
5. Further reading the road-map given by Mr.Ralf, I would like to develop the Bluestein's FFT algorithm.
https://gist.github.com/endolith/2783807 Regards Stéfan
![](https://secure.gravatar.com/avatar/5f88830d19f9c83e2ddfd913496c5025.jpg?s=120&d=mm&r=g)
On Fri, Jan 31, 2014 at 11:34 AM, Stéfan van der Walt <stefan@sun.ac.za>wrote:
Finding a suitable mentor for whatever project Jennifer chooses is an important factor in the choice of project, so I have to ask: do you have the bandwidth to be a mentor or help out this summer? Ralf
![](https://secure.gravatar.com/avatar/bae0c4adb1469b38cfe4b1844180c446.jpg?s=120&d=mm&r=g)
Thanks a ton Ralf, for bringing this up. It would my dream come true experience to work under any of you guys.However tough be the project, the guidance of a willing mentor shall make it a smooth and fruitful experience. Please let me know if you can help out. I would be really grateful to you. Regards
![](https://secure.gravatar.com/avatar/5f88830d19f9c83e2ddfd913496c5025.jpg?s=120&d=mm&r=g)
On Thu, Feb 6, 2014 at 11:59 AM, Stéfan van der Walt <stefan@sun.ac.za>wrote:
Great!
(of course, Ralf would be on top of that list :).
That depends a bit on the topic - I don't know much about scipy.special, Pauli is our resident guru.
Members of the dipy team would also be interested.
That's specifically for the spherical harmonics topic right? Ralf
![](https://secure.gravatar.com/avatar/af6c39d6943bd4b0e1fde23161e7bb8c.jpg?s=120&d=mm&r=g)
On 8 Feb 2014 04:51, "Ralf Gommers" <ralf.gommers@gmail.com> wrote:
Members of the dipy team would also be interested.
That's specifically for the spherical harmonics topic right?
Right. Spherical harmonics are used as bases in many of DiPy's reconstruction algorithms. You are right, though, that gsoc would also require an expert in special functions. Stéfan
![](https://secure.gravatar.com/avatar/5bc3911d30ec66d8a65a05846502f0c4.jpg?s=120&d=mm&r=g)
Hi Jennifer, On 31 January 2014 00:01, jennifer stone <jenny.stone125@gmail.com> wrote:
I once ported a method of Abate and Whitt to python. My aim was not to produce the nicest python implementation, but to stick closely to the code of Abate and Whitt in their paper. However, it might a useful starting point. http://nicky.vanforeest.com/queueing/euler/euler.html?highlight=laplace Nicky
![](https://secure.gravatar.com/avatar/96dd777e397ab128fedab46af97a3a4a.jpg?s=120&d=mm&r=g)
On Thu, Jan 30, 2014 at 4:01 PM, jennifer stone <jenny.stone125@gmail.com>wrote:
No comment, as I don't use this functionality. I don't know how many folks would want this.
This sounds very doable. How much work do you think would be involved?
4. Lastly, we are yet to have Inverse Laplace transforms which as Ralf has rightly pointed out it may be too challenging to implement.
This is more ambitious, I'm not in a position to comment on whether it is doable in the summer time frame.
5. Further reading the road-map given by Mr.Ralf, I would like to develop the Bluestein's FFT algorithm.
This one could be quite involved, but useful. The problem is not so much *a* Bluestein FFT, but combining it with the current FFTPACK so that factors other than 2,3,4, or 5 are handled with the Bluestein algorithm. FFTPACK is in Fortran and not very well documented. I wouldn't recommend this project unless you are pretty familiar with FFTs and Fortran. It is unfortunate that the latest versions of FFTPACK are GPL. A BSD licensed package that already implements the Bluestein algorithm for FFTs is Parallel Colt<https://sites.google.com/site/piotrwendykier/software/parallelcolt>, which is in Java but could maybe be translated. A similar but smaller project, not involving integration with the general FFT, would be a stand alone chirpz transform, might be too easy though ;)
Thanks for reading along till the end. I shall append to this mail as when I am struck with ideas. Please do give your valuable guidance
Chuck
![](https://secure.gravatar.com/avatar/bae0c4adb1469b38cfe4b1844180c446.jpg?s=120&d=mm&r=g)
3. As stated earlier, we have spherical harmonic functions (with much scope the second kind. I am confident about about the implementation of ellipsoidal H function of first kind but don't know much about the second kind. But I believe we can work it out in due course.And cylindrical harmonics can be carried out using Bessel functions.
![](https://secure.gravatar.com/avatar/da3a0a1942fbdc5ee9a9b8115ac5dae7.jpg?s=120&d=mm&r=g)
-----BEGIN PGP SIGNED MESSAGE----- Hash: SHA1 Hi, 04.02.2014 20:30, jennifer stone kirjoitti:
It's not so often someone wants to work on scipy.special, so you'd be welcome to improve it :) The general structure of work on special functions goes as follows: - - Check if there is a license-compatible implementation that someone has already written. This is usually not the case. - - Find formulas for evaluating the function in terms of more primitive operations. (Ie. power series, asymptotic series, continued fractions, expansions in terms of other special functions, ...) - - Determine the parameter region where the expansions converge in a floating point implementation, and select algorithms appropriately. Here it helps if you find a research paper where the author has already thought about what sort of an approach works best. - - Life is usually made *much* easier thanks to Fredrik Johansson's prior work on arbitrary-precision arithmetic library mpmath http://code.google.com/p/mpmath/ It can usually be used to check the "true" values of the functions. Also it contains implementations of algorithms for evaluating special functions, but because mpmath works with arbitrary precision numbers, these algorithms are not directly usable for floating-point calculations, as in floating point you cannot adjust the precision of the calculation dynamically. Moreover, the arbitrary-precision arithmetic can be slow compared to a more optimized floating point implementations. Spherical harmonics might be a reasonable part of a GSoC proposal. However, note that there exists also a *second* Legendre polynomial function `lpmv`, which doesn't store the values of the previous N functions. There's one numerical problem in the current way of evaluation via ~Pmn(cos(theta)), which is that this approach seems to lose relatively much precision at large orders for certain values of theta. I don't recall now exactly how imprecise it becomes at large orders, but it may be necessary to check. Adding new special functions also sounds like an useful project. Here, it helps if they are something that you expect you will need later on :) There's also the case that several of the functions in Scipy have only implementations for real-valued inputs, although the functions would be defined on the whole complex plane. A list of the situation is here: https://github.com/scipy/scipy/blob/master/scipy/special /generate_ufuncs.py#L85 Lowercase d correspond to real-valued implementations, uppercase D to complex-valued. I'm not at the moment completely sure which would have the highest priority --- whether you need this or not really depends on the application. If you want additional ideas about possible things to fix in scipy.special, take a look at this file: https://github.com/scipy/scipy/blob/master/scipy/special/tests /test_mpmath.py#L648 The entries marked @knownfailure* have some undiagnosed issues in the implementation, which might be useful to look into. However: most of these have to do with corner cases in hypergeometric functions. Trying to address those is likely a risky GSoC topic, as the multi-argument hyp* functions are challenging to evaluate in floating point. (mpmath and Mathematica can evaluate them in most parameter regimes, but AFAIK both require arbitrary-precision methods for this.) So I think there would be a large number of possible things to do here, and help would be appreciated. - -- Pauli Virtanen -----BEGIN PGP SIGNATURE----- Version: GnuPG v1.4.14 (GNU/Linux) iEYEARECAAYFAlL6iwAACgkQ6BQxb7O0pWBfOgCfYHAB12N4FWDmrqx8/ORTBRps pXYAoL3ufAiShe+0qTEGfEvrmDgr1X0p =kAwF -----END PGP SIGNATURE-----
![](https://secure.gravatar.com/avatar/bae0c4adb1469b38cfe4b1844180c446.jpg?s=120&d=mm&r=g)
On Wed, Feb 12, 2014 at 2:11 AM, Pauli Virtanen <pav@iki.fi> wrote:
values small N's by using recursion. Why can't we first check if m and n are n positive integers and pass them to lpmv itself rather than lpmn? lpmn does give us the derivatives too, but sph_harm has no need for that, and of course all the values for <m and <n too are trimmed. What is the benefit of lpmn over lpmv at present?
Adding new special functions also sounds like an useful project. Here, it helps if they are something that you expect you will need later on :)
I am unable to think beyond ellipsoidal functions. As for their use, we as students used them in problems of thermal equilibrium in ellipsoidal bodies, and some scattering cases. Though I have no idea if its use in general is quite prominent. And cylindrical harmonic function would be useful but I feel it's quite straight forward to implement (correct me if I am wrong) .
Yeah, many of the known failures seem to revolve around hyp2f1. An unexplained inclination towards hypergeometric functions really tempts me to plunge into this. If it's too risky, I can work on this after the summers, as I would have gained quite a lot of experience with the code here. So I think there would be a large number of possible things to do
Regard Jennifer
![](https://secure.gravatar.com/avatar/da3a0a1942fbdc5ee9a9b8115ac5dae7.jpg?s=120&d=mm&r=g)
16.02.2014 23:34, Jennifer stone kirjoitti: [clip]
If you are interested in the hypergeometric numerical evaluation, it's probably a good idea to take a look at this recent master's thesis written on the problem: http://people.maths.ox.ac.uk/porterm/research/pearson_final.pdf This may give some systematic overview on the range of methods available. (Note that for copyright reasons, it's not a good idea to look closely at the source codes linked from that thesis, as they are not available under a compatible license.) It may well be that the best approach for evaluating these functions, if accuracy in the whole parameter range is wanted, in the end turns out to require arbitrary-precision computations. In that case, it would be a very good idea to look at how the problem is approached in mpmath. There are existing multiprecision packages written in C, and using one of them in scipy.special could bring better evaluation performance even if the algorithm is the same. -- Pauli Virtanen
![](https://secure.gravatar.com/avatar/bae0c4adb1469b38cfe4b1844180c446.jpg?s=120&d=mm&r=g)
If you are interested in the hypergeometric numerical evaluation, it's
Yeah, this seems to be brilliant idea. mpmath too, I assume, must have used some of the methods mentioned in the thesis. I ll look through the code and get back. I am still unaware of the complexity of project expected at GSoC. This project looks engaging to me. Will an attempt to improve both Spherical harmonic functions ( improving the present algorithm to avoid the calculation for lower n's and m's) and hypergeometric functions be too ambitious or is it doable? Regards Jennifer
-- Pauli Virtanen
![](https://secure.gravatar.com/avatar/bae0c4adb1469b38cfe4b1844180c446.jpg?s=120&d=mm&r=g)
In an attempt to analyze the accuracy of hyp2f1, Different cases mentioned in Abramowitz ( http://people.math.sfu.ca/~cbm/aands/page_561.htm<http://bl-1.com/click/load/UmcPPANhUGxVMQNuUGU-b0231> ) and also in the Thesis on 'Computation of Hypergeometric functions" (http://people.maths.ox.ac.uk/porterm/research/pearson_final.pdf<http://bl-1.com/click/load/VmMBMlw-b0221ADxRNVI-b0169BDY-b0231>, pg 65-66) were tried out, and the function fails without warning when: c<0, c is not integral |c|>>|a| and |b| For example: sp.hyp2f1(10,5,-300.5,0.5)
-6.5184949735e+156 while the answer is
*-3.8520770815e+32* this case appears to filter down to hys2f1 in the source code (scipy.special.cephes.hyp2f1) I tried the same input in mpmath to check if it works there: hyp2f1(10,5,-300.5,0.5)
mpf('0.9211827166328477893913199888') which is the solution when we apply power series expansion.
however MATLAB succeeds in giving the required solution. Another interesting fact is that of the methods mentioned in the thesis: Taylor series expansion, fraction method with double precision, Gauss-Jacobi method and RK4), none succeeds in the given case. I don't have any idea how the function itself is evaluated in the given case. Any leads on how it is done and how MATLAB executes it? On Thu, Feb 20, 2014 at 1:16 AM, Jennifer stone <jenny.stone125@gmail.com>wrote:
![](https://secure.gravatar.com/avatar/da3a0a1942fbdc5ee9a9b8115ac5dae7.jpg?s=120&d=mm&r=g)
23.02.2014 15:18, Jennifer stone kirjoitti: [clip]
Typically, you need to tell mpmath to use an appropriate precision for the evaluation:
You can look in their documentation what references the implementation is based on. However, it's not a good idea to go beyond that in wondering how it does things --- from legal POV and from fair play. *** What can help in hyp2f1 for large values of a,b,c is use of recurrence relations. These are typically stable in one direction only. [1] This seems to be still a partially open research question... Our current hyp2f1 implementation does use recurrences (hyp2f1ra), but perhaps they are not invoked for this case. The problem here can be the accurate determination of the convergence region for each parameter value. [1] http://www.ams.org/journals/mcom/2007-76-259/S0025-5718-07-01918-7/ -- Pauli Virtanen
![](https://secure.gravatar.com/avatar/bae0c4adb1469b38cfe4b1844180c446.jpg?s=120&d=mm&r=g)
sp.hyp2f1(10,5,-300.5,0.5)
What are the implications of adopting this method in the cases where hyp2f1 fails to give accurate results? Except of course the fact that this implementation would be heavy. direct implementation without a change in sensitivity fails to give the required answer. Regards Jenny
-- Pauli Virtanen
![](https://secure.gravatar.com/avatar/5f88830d19f9c83e2ddfd913496c5025.jpg?s=120&d=mm&r=g)
On Fri, Jan 31, 2014 at 12:01 AM, jennifer stone <jenny.stone125@gmail.com>wrote:
Another idea: add support for discrete wavelet transforms in scipy.signal. There's a fair bit of interest for those here I think. It would start by integrating https://github.com/rgommers/pywt, then adding some new features. Feature ideas: - 1-D and 2-D inverse SWT (have been requested several times on the PyWavelets list and issue tracker). - signal denoising (SureShrink & co, for scipy.signal) - image compression/denoising/... algorithms (for scikit-image) Ralf
![](https://secure.gravatar.com/avatar/af6c39d6943bd4b0e1fde23161e7bb8c.jpg?s=120&d=mm&r=g)
On Fri, 31 Jan 2014 04:31:01 +0530, jennifer stone wrote:
As stated before, I am personally interested in seeing the spherical harmonics in SciPy improve.
5. Further reading the road-map given by Mr.Ralf, I would like to develop the Bluestein's FFT algorithm.
https://gist.github.com/endolith/2783807 Regards Stéfan
![](https://secure.gravatar.com/avatar/5f88830d19f9c83e2ddfd913496c5025.jpg?s=120&d=mm&r=g)
On Fri, Jan 31, 2014 at 11:34 AM, Stéfan van der Walt <stefan@sun.ac.za>wrote:
Finding a suitable mentor for whatever project Jennifer chooses is an important factor in the choice of project, so I have to ask: do you have the bandwidth to be a mentor or help out this summer? Ralf
![](https://secure.gravatar.com/avatar/bae0c4adb1469b38cfe4b1844180c446.jpg?s=120&d=mm&r=g)
Thanks a ton Ralf, for bringing this up. It would my dream come true experience to work under any of you guys.However tough be the project, the guidance of a willing mentor shall make it a smooth and fruitful experience. Please let me know if you can help out. I would be really grateful to you. Regards
![](https://secure.gravatar.com/avatar/5f88830d19f9c83e2ddfd913496c5025.jpg?s=120&d=mm&r=g)
On Thu, Feb 6, 2014 at 11:59 AM, Stéfan van der Walt <stefan@sun.ac.za>wrote:
Great!
(of course, Ralf would be on top of that list :).
That depends a bit on the topic - I don't know much about scipy.special, Pauli is our resident guru.
Members of the dipy team would also be interested.
That's specifically for the spherical harmonics topic right? Ralf
![](https://secure.gravatar.com/avatar/af6c39d6943bd4b0e1fde23161e7bb8c.jpg?s=120&d=mm&r=g)
On 8 Feb 2014 04:51, "Ralf Gommers" <ralf.gommers@gmail.com> wrote:
Members of the dipy team would also be interested.
That's specifically for the spherical harmonics topic right?
Right. Spherical harmonics are used as bases in many of DiPy's reconstruction algorithms. You are right, though, that gsoc would also require an expert in special functions. Stéfan
![](https://secure.gravatar.com/avatar/5bc3911d30ec66d8a65a05846502f0c4.jpg?s=120&d=mm&r=g)
Hi Jennifer, On 31 January 2014 00:01, jennifer stone <jenny.stone125@gmail.com> wrote:
I once ported a method of Abate and Whitt to python. My aim was not to produce the nicest python implementation, but to stick closely to the code of Abate and Whitt in their paper. However, it might a useful starting point. http://nicky.vanforeest.com/queueing/euler/euler.html?highlight=laplace Nicky
![](https://secure.gravatar.com/avatar/96dd777e397ab128fedab46af97a3a4a.jpg?s=120&d=mm&r=g)
On Thu, Jan 30, 2014 at 4:01 PM, jennifer stone <jenny.stone125@gmail.com>wrote:
No comment, as I don't use this functionality. I don't know how many folks would want this.
This sounds very doable. How much work do you think would be involved?
4. Lastly, we are yet to have Inverse Laplace transforms which as Ralf has rightly pointed out it may be too challenging to implement.
This is more ambitious, I'm not in a position to comment on whether it is doable in the summer time frame.
5. Further reading the road-map given by Mr.Ralf, I would like to develop the Bluestein's FFT algorithm.
This one could be quite involved, but useful. The problem is not so much *a* Bluestein FFT, but combining it with the current FFTPACK so that factors other than 2,3,4, or 5 are handled with the Bluestein algorithm. FFTPACK is in Fortran and not very well documented. I wouldn't recommend this project unless you are pretty familiar with FFTs and Fortran. It is unfortunate that the latest versions of FFTPACK are GPL. A BSD licensed package that already implements the Bluestein algorithm for FFTs is Parallel Colt<https://sites.google.com/site/piotrwendykier/software/parallelcolt>, which is in Java but could maybe be translated. A similar but smaller project, not involving integration with the general FFT, would be a stand alone chirpz transform, might be too easy though ;)
Thanks for reading along till the end. I shall append to this mail as when I am struck with ideas. Please do give your valuable guidance
Chuck
![](https://secure.gravatar.com/avatar/bae0c4adb1469b38cfe4b1844180c446.jpg?s=120&d=mm&r=g)
3. As stated earlier, we have spherical harmonic functions (with much scope the second kind. I am confident about about the implementation of ellipsoidal H function of first kind but don't know much about the second kind. But I believe we can work it out in due course.And cylindrical harmonics can be carried out using Bessel functions.
![](https://secure.gravatar.com/avatar/da3a0a1942fbdc5ee9a9b8115ac5dae7.jpg?s=120&d=mm&r=g)
-----BEGIN PGP SIGNED MESSAGE----- Hash: SHA1 Hi, 04.02.2014 20:30, jennifer stone kirjoitti:
It's not so often someone wants to work on scipy.special, so you'd be welcome to improve it :) The general structure of work on special functions goes as follows: - - Check if there is a license-compatible implementation that someone has already written. This is usually not the case. - - Find formulas for evaluating the function in terms of more primitive operations. (Ie. power series, asymptotic series, continued fractions, expansions in terms of other special functions, ...) - - Determine the parameter region where the expansions converge in a floating point implementation, and select algorithms appropriately. Here it helps if you find a research paper where the author has already thought about what sort of an approach works best. - - Life is usually made *much* easier thanks to Fredrik Johansson's prior work on arbitrary-precision arithmetic library mpmath http://code.google.com/p/mpmath/ It can usually be used to check the "true" values of the functions. Also it contains implementations of algorithms for evaluating special functions, but because mpmath works with arbitrary precision numbers, these algorithms are not directly usable for floating-point calculations, as in floating point you cannot adjust the precision of the calculation dynamically. Moreover, the arbitrary-precision arithmetic can be slow compared to a more optimized floating point implementations. Spherical harmonics might be a reasonable part of a GSoC proposal. However, note that there exists also a *second* Legendre polynomial function `lpmv`, which doesn't store the values of the previous N functions. There's one numerical problem in the current way of evaluation via ~Pmn(cos(theta)), which is that this approach seems to lose relatively much precision at large orders for certain values of theta. I don't recall now exactly how imprecise it becomes at large orders, but it may be necessary to check. Adding new special functions also sounds like an useful project. Here, it helps if they are something that you expect you will need later on :) There's also the case that several of the functions in Scipy have only implementations for real-valued inputs, although the functions would be defined on the whole complex plane. A list of the situation is here: https://github.com/scipy/scipy/blob/master/scipy/special /generate_ufuncs.py#L85 Lowercase d correspond to real-valued implementations, uppercase D to complex-valued. I'm not at the moment completely sure which would have the highest priority --- whether you need this or not really depends on the application. If you want additional ideas about possible things to fix in scipy.special, take a look at this file: https://github.com/scipy/scipy/blob/master/scipy/special/tests /test_mpmath.py#L648 The entries marked @knownfailure* have some undiagnosed issues in the implementation, which might be useful to look into. However: most of these have to do with corner cases in hypergeometric functions. Trying to address those is likely a risky GSoC topic, as the multi-argument hyp* functions are challenging to evaluate in floating point. (mpmath and Mathematica can evaluate them in most parameter regimes, but AFAIK both require arbitrary-precision methods for this.) So I think there would be a large number of possible things to do here, and help would be appreciated. - -- Pauli Virtanen -----BEGIN PGP SIGNATURE----- Version: GnuPG v1.4.14 (GNU/Linux) iEYEARECAAYFAlL6iwAACgkQ6BQxb7O0pWBfOgCfYHAB12N4FWDmrqx8/ORTBRps pXYAoL3ufAiShe+0qTEGfEvrmDgr1X0p =kAwF -----END PGP SIGNATURE-----
![](https://secure.gravatar.com/avatar/bae0c4adb1469b38cfe4b1844180c446.jpg?s=120&d=mm&r=g)
On Wed, Feb 12, 2014 at 2:11 AM, Pauli Virtanen <pav@iki.fi> wrote:
values small N's by using recursion. Why can't we first check if m and n are n positive integers and pass them to lpmv itself rather than lpmn? lpmn does give us the derivatives too, but sph_harm has no need for that, and of course all the values for <m and <n too are trimmed. What is the benefit of lpmn over lpmv at present?
Adding new special functions also sounds like an useful project. Here, it helps if they are something that you expect you will need later on :)
I am unable to think beyond ellipsoidal functions. As for their use, we as students used them in problems of thermal equilibrium in ellipsoidal bodies, and some scattering cases. Though I have no idea if its use in general is quite prominent. And cylindrical harmonic function would be useful but I feel it's quite straight forward to implement (correct me if I am wrong) .
Yeah, many of the known failures seem to revolve around hyp2f1. An unexplained inclination towards hypergeometric functions really tempts me to plunge into this. If it's too risky, I can work on this after the summers, as I would have gained quite a lot of experience with the code here. So I think there would be a large number of possible things to do
Regard Jennifer
![](https://secure.gravatar.com/avatar/da3a0a1942fbdc5ee9a9b8115ac5dae7.jpg?s=120&d=mm&r=g)
16.02.2014 23:34, Jennifer stone kirjoitti: [clip]
If you are interested in the hypergeometric numerical evaluation, it's probably a good idea to take a look at this recent master's thesis written on the problem: http://people.maths.ox.ac.uk/porterm/research/pearson_final.pdf This may give some systematic overview on the range of methods available. (Note that for copyright reasons, it's not a good idea to look closely at the source codes linked from that thesis, as they are not available under a compatible license.) It may well be that the best approach for evaluating these functions, if accuracy in the whole parameter range is wanted, in the end turns out to require arbitrary-precision computations. In that case, it would be a very good idea to look at how the problem is approached in mpmath. There are existing multiprecision packages written in C, and using one of them in scipy.special could bring better evaluation performance even if the algorithm is the same. -- Pauli Virtanen
![](https://secure.gravatar.com/avatar/bae0c4adb1469b38cfe4b1844180c446.jpg?s=120&d=mm&r=g)
If you are interested in the hypergeometric numerical evaluation, it's
Yeah, this seems to be brilliant idea. mpmath too, I assume, must have used some of the methods mentioned in the thesis. I ll look through the code and get back. I am still unaware of the complexity of project expected at GSoC. This project looks engaging to me. Will an attempt to improve both Spherical harmonic functions ( improving the present algorithm to avoid the calculation for lower n's and m's) and hypergeometric functions be too ambitious or is it doable? Regards Jennifer
-- Pauli Virtanen
![](https://secure.gravatar.com/avatar/bae0c4adb1469b38cfe4b1844180c446.jpg?s=120&d=mm&r=g)
In an attempt to analyze the accuracy of hyp2f1, Different cases mentioned in Abramowitz ( http://people.math.sfu.ca/~cbm/aands/page_561.htm<http://bl-1.com/click/load/UmcPPANhUGxVMQNuUGU-b0231> ) and also in the Thesis on 'Computation of Hypergeometric functions" (http://people.maths.ox.ac.uk/porterm/research/pearson_final.pdf<http://bl-1.com/click/load/VmMBMlw-b0221ADxRNVI-b0169BDY-b0231>, pg 65-66) were tried out, and the function fails without warning when: c<0, c is not integral |c|>>|a| and |b| For example: sp.hyp2f1(10,5,-300.5,0.5)
-6.5184949735e+156 while the answer is
*-3.8520770815e+32* this case appears to filter down to hys2f1 in the source code (scipy.special.cephes.hyp2f1) I tried the same input in mpmath to check if it works there: hyp2f1(10,5,-300.5,0.5)
mpf('0.9211827166328477893913199888') which is the solution when we apply power series expansion.
however MATLAB succeeds in giving the required solution. Another interesting fact is that of the methods mentioned in the thesis: Taylor series expansion, fraction method with double precision, Gauss-Jacobi method and RK4), none succeeds in the given case. I don't have any idea how the function itself is evaluated in the given case. Any leads on how it is done and how MATLAB executes it? On Thu, Feb 20, 2014 at 1:16 AM, Jennifer stone <jenny.stone125@gmail.com>wrote:
![](https://secure.gravatar.com/avatar/da3a0a1942fbdc5ee9a9b8115ac5dae7.jpg?s=120&d=mm&r=g)
23.02.2014 15:18, Jennifer stone kirjoitti: [clip]
Typically, you need to tell mpmath to use an appropriate precision for the evaluation:
You can look in their documentation what references the implementation is based on. However, it's not a good idea to go beyond that in wondering how it does things --- from legal POV and from fair play. *** What can help in hyp2f1 for large values of a,b,c is use of recurrence relations. These are typically stable in one direction only. [1] This seems to be still a partially open research question... Our current hyp2f1 implementation does use recurrences (hyp2f1ra), but perhaps they are not invoked for this case. The problem here can be the accurate determination of the convergence region for each parameter value. [1] http://www.ams.org/journals/mcom/2007-76-259/S0025-5718-07-01918-7/ -- Pauli Virtanen
![](https://secure.gravatar.com/avatar/bae0c4adb1469b38cfe4b1844180c446.jpg?s=120&d=mm&r=g)
sp.hyp2f1(10,5,-300.5,0.5)
What are the implications of adopting this method in the cases where hyp2f1 fails to give accurate results? Except of course the fact that this implementation would be heavy. direct implementation without a change in sensitivity fails to give the required answer. Regards Jenny
-- Pauli Virtanen
![](https://secure.gravatar.com/avatar/5f88830d19f9c83e2ddfd913496c5025.jpg?s=120&d=mm&r=g)
On Fri, Jan 31, 2014 at 12:01 AM, jennifer stone <jenny.stone125@gmail.com>wrote:
Another idea: add support for discrete wavelet transforms in scipy.signal. There's a fair bit of interest for those here I think. It would start by integrating https://github.com/rgommers/pywt, then adding some new features. Feature ideas: - 1-D and 2-D inverse SWT (have been requested several times on the PyWavelets list and issue tracker). - signal denoising (SureShrink & co, for scipy.signal) - image compression/denoising/... algorithms (for scikit-image) Ralf
participants (7)
-
Charles R Harris
-
jennifer stone
-
Jennifer stone
-
nicky van foreest
-
Pauli Virtanen
-
Ralf Gommers
-
Stéfan van der Walt