[SciPy-dev] Volunteer for Scipy Project
Charles R Harris
Mon Oct 5 23:36:10 CDT 2009
On Mon, Oct 5, 2009 at 10:14 PM, Anne Archibald
> 2009/10/5 Charles R Harris <firstname.lastname@example.org>:
> > On Mon, Oct 5, 2009 at 7:29 PM, Anne Archibald <
> > wrote:
> >> 2009/10/5 Charles R Harris <email@example.com>:
> >> >
> >> >
> >> > On Mon, Oct 5, 2009 at 3:20 PM, Anne Archibald
> >> > <firstname.lastname@example.org>
> >> > wrote:
> >> >>
> >> >> 2009/10/5 David Goldsmith <email@example.com>:
> >> >> > Just curious, Anne: have you anything in particular in mind (i.e.,
> >> >> > are
> >> >> > there
> >> >> > some small - or gaping - holes in scipy (IYO, of course) which you
> >> >> > know
> >> >> > could be filled by a careful implementation of something(s) extant
> >> >> > the
> >> >> > literature)?
> >> >>
> >> >> Well, not exactly - the examples I had in mind were minor and/or in
> >> >> the past. I ran into problems with scipy's hyp2f1, for example, so I
> >> >> went and looked up the best algorithm I could find for it (and I
> >> >> I contributed that code). I wanted the Kuiper test as an alternative
> >> >> to the Kolmogorov-Smirnov test (it's invariant under cyclic
> >> >> permutations, and is sensitive to different features of the
> >> >> distribution) so I looked up the test and the special function needed
> >> >> to interpret its results. (I haven't contributed this to scipy yet,
> >> >> mostly because I chose an interface that's not particularly
> >> >> with that for scipy's K-S test.) And on a larger scale, that's what
> >> >> scipy.spatial's kdtree implementation is.
> >> >>
> >> >> For examples where I think a bit of lit review plus implementation
> >> >> work might help, I'd say that the orthogonal polynomials could use
> >> >> some work - the generic implementation in scipy.special falls apart
> >> >> rapidly as you go to higher orders. I always implement my own
> >> >> Chebyshev polynomials using the cos(n*arccos(x)) expression, for
> >> >> example, and special implementations for the others might be very
> >> >> useful.
> >> >>
> >> >
> >> > At what order does the scipy implementation of the Chebyshev
> >> > fall apart? I looked briefly at that package a long time ago, but
> >> > used
> >> > it. I ask so I can check the chebyshev module that is going into
> >> By n=30 they are off by as much as 0.0018 on [-1,1]; n=31 they are off
> >> by 0.1, and by n=35 they are off by four - not great for values that
> >> should be in the interval [-1,1]. This may seem like an outrageously
> >> high degree for a polynomial, but there's no reason they need be this
> >> bad, and one could quite reasonably want to use an order this high,
> >> say for function approximation.
> > Hmm, I get an maximum error of about 1e-13 for n=100 using my routine,
> > isn't great and can be improved a bit with a few tricks, but is probably
> > good enough. That's using simple Clenshaw recursion for the evaluation.
> > There must be something seriously wrong with the scipy version. I
> > use fits with n > 50 because truncating such a series gives a good
> > approximation to the minmax fit and it's also nice to see how the
> > coefficients fall off to estimate the size of the truncation error.
> Upon closer inspection, it seems that scipy starts from the recurrence
> relation but then uses it only to extract roots, which it then uses to
> construct a poly1d instance. (Which presumably then goes on to convert
> it into the 1, x, x**2, ... basis.) This is numerically disastrous.
> The reason it does this is because scipy doesn't (didn't) have a
> general-purpose polynomial class that can keep the polynomial in other
> representations. What do you think of the idea of a polynomial class
> that supports several different bases for the set of polynomials? (At
> least, the power basis and Chebyshev polynomials, but possibly also
> the other orthogonal polynomials. Translated and scaled bases might
> come in handy too.)
The new Chebyshev class has a domain [beg, end] that defines the translation
and scaling. The polynomial class needs one also. In fact, I would like to
define a Polynomial class to replace poly1d.
> Too general? Not general enough (barycentric
> representations, representation in terms of roots)? I think it could
> be done without too much difficulty.
The difficulty is in between. I've been thinking of a base series class that
would handle the domains and access to the coefficients, as well as
addition/subtraction, scalar multiplication/division, and pickling. Derived
classes could then overload the
division/multiplication/integration/differentiation. Roots would come
directly from the recursion relation. I don't know if the stability of the
Clenshaw recursion (essentially synthetic division with the value as the
remainder) is guaranteed for the polynomials/domains of interest, but we
could make a list and look into it.
-------------- next part --------------
An HTML attachment was scrubbed...
More information about the Scipy-dev