[Numpy-discussion] Suggestions for GSoC Projects

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

Hash: SHA1


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

  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


  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:


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:


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
Version: GnuPG v1.4.14 (GNU/Linux)


More information about the NumPy-Discussion mailing list