[SciPy-user] Polynomial interpolation

[Rob Clewley]

On Mon, Apr 21, 2008 at 6:07 AM, Joris De Ridder
<Joris.DeRidder at ster.kuleuven.be> wrote:
>
>  On 21 Apr 2008, at 08:02, Stéfan van der Walt wrote:
>
>  > Talking about derivatives, does anyone know whether
>  >
>  > http://en.wikipedia.org/wiki/Automatic_differentiation
>  >
>  > is of value?  It's been on my TODO list for a while, but I haven't
>  > gotten round to studying it in detail.

It is of great value for several reasons (e.g. see Eric Phipps' thesis
at www.math.cornell.edu/~gucken/PDF/phipps.pdf). Our group at Cornell
a couple of years ago had been waiting to see if there would be a
standard package emerging for us to interface into Python....

>  I think the Scientific.Functions.Derivatives subpackage  of Konrad
>  Hinsen has such functionality.

For one thing, I strongly dislike the way of interacting with the
numeric objects and their derivatives in that package. It's not very
Pythonic in its use of classes and duck typing.

The other problem is, the situations where it would be of most use
need so many computations that it's not an effective approach unless
it is done at the level of C code. A couple of years ago we tried
interfacing to a new release of ADOL-C through SWIG, but found some
extremely strange memory errors that we couldn't sort out. I think
those only showed up when we tried to pass in data that had been
prepared by SWIG from Python arrays. Anyway, getting a package like
that properly interfaced would be the way forward as far as we are
concerned. OpenAD looks like another good bet for an open-source
library. Then again, the problem of efficiency becomes getting the
users functions into C code so that "source code transformation" can
be performed.

I certainly like the idea of having all in-built functions "knowing"
their derivatives, but it's not clear how these python-level
representations can be best interfaced to C code, whether the basis
for the AD is "source code transformation" or "operator overloading".
I think there would need to be a new class that allows "user"
functions that know their derivatives but which are defined in a
piecewise-fashion, e.g. to include solutions to differential equations
(for instance) represented as interpolated polynomials.

>  Also, to come back to the original
>  thread subject,  Scientific.Functions.Interpolation is also able to
>  generate Interpolators (don't know the algorithm).

It's only linear interpolation, but, on the plus side, does support
multi-dimensional meshes, which is a generalization I wholeheartedly
endorse. Alas, multivariate polynomials or wavelet-type bases would be
needed.

<soapbox>
If we're going to start thinking "big" for supporting more
mathematically natural functionality, I believe we ought to be
thinking far enough out to support the basic objects of PDE
computations too (or at least a compatible class structure and API),
even if it's not fully utilized just yet. Scipy should support
scientific computation based around mathematically-oriented
fundamental objects and functionality (i.e. to hide the "dumbness" of
arrays inside some sugar coating). I think writing and debugging
complex mathematical tools will be a lot easier if we raise the bar a
little and use somewhat more sophisticated basic objects than arrays
(our efforts with Pointsets have helped enormously in writing
sub-packages of PyDSTool, in particular for putting PyCont together so
quickly). Such an approach will also be a lot easier to introduce to
new and less technical scientific users. The field of scientific
computation is still weighed down by keeping old Fortran-style APIs up
to date for re-use (of course, I'm guilty of this). Python brings a
fresh opportunity to break these shackles, at least interpreted in
terms of adding an extra level of indirection and duck-typed magic at
the heart of scipy. Some imagination is needed here (a la "import
antigravity" ;).
</soapbox>

-Rob