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

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

_______________________________________________

NumPy-Discussion mailing list

[hidden email]
http://mail.scipy.org/mailman/listinfo/numpy-discussion