[SciPy-dev] Generic polynomials class (was Re: Volunteer for Scipy Project)

[Anne Archibald]

2009/10/15 Fernando Perez <fperez.net at gmail.com>:
> On Thu, Oct 15, 2009 at 6:36 PM, Charles R Harris
> <charlesr.harris at gmail.com> wrote:
>>
>> How did you build the analysis objects? Part of what we have been discussing
>> it the easiest way to expose a consistent interface while reusing the
>> maximum amount of code.
>
> See
>
> http://github.com/fperez/nitime/blob/master/nitime/timeseries.py
>
> starting around line 504.  The analyzers store a reference to a time
> series object and then expose properties (that are replaced in-place
> by a static attribute once computed, as explained in the paper) that
> make the library calls.  All the actual numerics are in
>
> http://github.com/fperez/nitime/blob/master/nitime/algorithms.py
>
> I'm not saying this will exactly solve your current design issues, but
> it may provide some useful ideas to think about.

Well, I think I'm convinced of the necessity for a low-level
procedural interface, so I've set up a draft one in git. There's some
question of just how bare-bones it should really be - for example
should it include an add function, since that's just using the
coefficient-extension function then adding? What about multiplying,
which is the same thing, though this is not obvious?

As it stands I've only included implementations of functions with
nontrivial algorithms. This has the advantage that, realistically,
these functions will not be shared between different polynomial
representations. If I start including more straight-forward functions,
I will need to figure out how to share the implementation of "add"
between different representations if I want to avoid it being
boilerplate code for all representations.

If we're implementing the heavy lifting in a procedural interface,
then the object implementations will just be plumbing (as I see your
analyzer objects are). So on the one hand there's not much code to be
shared, but on the other it's all boilerplate that would benefit from
code sharing.

I think, though, that the next important decision has nothing direct
to do with code sharing; we need to settle on what polynomial objects
will look like to users.

Some things are clear:

* Polynomial objects will be immutable.
* Polynomial objects will support + - * / % ** divmod and __call__
* Polynomial objects will be able to report their roots
* Polynomial objects will be able to return polynomials representing
their derivatives and antiderivatives
* Some polynomial objects will have extra features, like the ability
to trim small coefficients from Chebyshev polynomials
* We need some method of converting from one representation to another
* We should have convenience objects zero, one, and X, to build polynomials.
* There should be a least-squares polynomial fitting function.
* Some polynomials will have functions that generate them (e.g.
Chebyshev series for a function)
* We need at least power basis, Chebyshev, and Lagrange polynomials.
* Ultimately the polynomial class should be able to replace
KroghInterpolator and BarycentricInterpolator and be stably returned
from the orthogonal polynomial routines

Less clear:

* Should polynomials store a degree attribute or should it be inferred
from the size of their coefficient array? (Sometimes it's handy to
work with degree-inflated polynomials, but deflating their degree is a
nontrivial operation.)
* Should we try to accommodate complex variables? Complex coefficients?
* Should we try to accommodate vector values?
* Should polynomials have a method to evaluate derivatives and
antiderivatives directly without constructing the
derivative/antiderivative object first?
* Do we need Lagrange-Hermite (some derivatives specified as well as
values) polynomials? (KroghInterpolator implements these)
* How should code that works for any polynomial type (e.g.
multiplication by a Lagrange polynomial) identify that an object it
has been given is some sort of polynomial? Callable and presence of a
degree attribute? Subclass of Polynomial?
* How should users implement their own polynomial types so they
interoperate with scipy's built-in ones?

Still debated:

* I think we need a first-class object to represent the basis used to
express a polynomial. Charles Harris disagrees; he wants to encode
that information in the class of the polynomial plus an interval
stored in the polynomial. Some different scheme would necessarily be
used for Lagrange polynomials.
* How should the zero polynomial be represented? I think its
coefficient array should be np.zeros(0); Charles Harris prefers
np.zeros(1).


Anne