# [Numpy-discussion] Suggestions for GSoC Projects

Pauli Virtanen pav at iki.fi
Tue Feb 11 15:41:41 EST 2014

```-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1

Hi,

04.02.2014 20:30, jennifer stone kirjoitti:
> 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.
>>
>> This sounds very doable. How much work do you think would be
>> involved?
>
> As Stefan so rightly pointed out, the function for spherical
> harmonic function, sph_harm at present calls lpmn thus evaluating
> all orders <N. An initial glance at the code and the algorithm
> gives me a feeling that it would be very well possible to avoid
> that by maybe avoiding the dependence on lpmn.
>
> Further, we can introduce ellipsoidal harmonic functions of first
> kind and 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.

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

- - Life is usually made *much* easier thanks to Fredrik Johansson's
prior work on arbitrary-precision arithmetic library 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.

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-----

```