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

[Anne Archibald]

2009/10/14 Charles R Harris <charlesr.harris at gmail.com>:
>
>
> On Tue, Oct 13, 2009 at 11:03 PM, Anne Archibald <peridot.faceted at gmail.com>
> wrote:
>>
>> 2009/10/13 Charles R Harris <charlesr.harris at gmail.com>:
>> >
>> > On Tue, Oct 13, 2009 at 7:24 PM, Anne Archibald
>> > <peridot.faceted at gmail.com>
>> > wrote:
>> >> specifying the interval. And if people do want to specify the
>> >> interval, using a basis object makes it much easier to avoid
>> >> forgetting it and getting nonsensical results (or worse, sensible
>> >> results until you try a different interval).
>> >>
>> >
>> > That is because they are low level building blocks, they should have the
>> > absolute minimum of crud stuck on, and that means the interval is [-1,
>> > 1].
>> > If you need more, that is what a class is for. That how low level things
>> > should work, small functions doing the minimum amount provide the most
>> > flexibility. The programmer should be handed enough rope to hang
>> > themselves
>> > at that level if they so choose. That's why folks use C and not Pascal.
>>
>> Well, the glib answer would be that we're using python, so we should
>> choose a pythonic way of doing things. But clearly having the
>> low-level functions is important to some potential users - you and
>> Robert Kern, plus surely more other people - so where possible the
>> code should accommodate the as-low-level-as-possible approach.
>>
>> How's this for a compromise? The high-level class interface has a
>> Polynomial object and various Basis subclasses. For representations
>> where having a featureful Basis object is important, the
>> implementation should just write one, inheriting as appropriate. For
>> representations that don't need this extra stuff, people can just
>> write a collection of functions acting on coefficient arrays (as you
>> did for the Chebyshev polynomials) and create a basis object using
>>
>
> I'm trying to get you to consider not using inheritance at all and see how
> things look. We might not want to go that way in the end, but inheritance
> tends to be overused. Probably because it is the first tool that comes to
> mind and because it looks like a cheap way to reuse code. But it is often
> not the best tool so folks need to step back and see if there aren't better
> ways to get things done. The fact that you are putting some functionality in
> the base class indicates that you are thinking along the lines of
> implementation. Passing a Basis object into a constructor also doesn't smell
> right. The question to ask  is what role you want the Basis object to play
> and if it doesn't more naturally belong somewhere else if you weren't
> constrained by inheritance. In a sense it looks like you would like to
> derive the interface class from the basis object, but somehow that doesn't
> seem to work well with inheritance because then it is hard for the derived
> objects to share code.

Well, conceptually there are two different objects for a given
representation: a basis, and a polynomial represented as a linear
combination of basis elements.

The polynomial object should obviously have methods for things like
multiplication and addition. I agree with your idea that it should be
immutable. But it must have some way of knowing what basis its
coefficient refers to.

The basis object is a natural place to put information like roots of
basis functions, weights for barycentric interpolation, and derivative
matrices. Having basis objects is also a natural way to answer the
question "are these polynomials represented in the same basis?",
necessary for type-checking; it also lets user code create new objects
in the same basis as a polynomial it is handed. Being able to define a
custom __eq__ allows users to create bases in different places that
are nonetheless compatible. (Incidentally, my intent is that users
should not ever call the Polynomial constructor; the preferred way to
make a polynomial object is to hand a basis object a coefficient
array.)

Given that we need these two kinds of object, how should the code be
arranged? We have, so far, three different representations of
polynomials to work with, so we can look at how much code can
reasonably be shared between them.

Addition can be identical for all three, provided that we have an
operation to "extend" a coefficient array to a basis with more
elements. This extension operation is identical - simply padding with
zeros - for Chebyshev and power bases, but different for the Lagrange
basis.

Evaluation is different for all three bases (but will be similar for
the various orthogonal-polynomial bases if we implement them).

Multiplication and division are different for all three bases, though
I know how to write a generic division that should both work and be
numerically stable (albeit not as efficient as one might wish) for all
representations.

Power can be implemented once for all types, but has a more efficient
specialized implementation for Lagrange polynomials.

Companion matrices and roots can be done generically for any
representation that allows division, but will have special
implementations in most bases.

All the various adapter plumbing to make operators work on polynomials
is identical for all representations.

The convert() operation to change the representation of a polynomial,
can probably work fine given the destination basis and a polynomial
object, but will need to be different for each basis.

So to summarize: we have some code that should be shared by all
representations, some that can be written generically but some
representations will provide specialized implementations, some that is
shared between only some representations, and some that is unique to
each representation. User representations, if they occur, may share
some code with existing representations, but it's not clear which.

So where do we put the code? Whether any given bit of code should be
in the polynomial objects or the basis objects is not clear; some bits
need to be in the basis objects (conversion code, zero(), one(), X())
and some bits need to be in the polynomial objects (type checking,
operator implementation). But the actual meat of a new representation
could in principle go in either place or in some third place (module
namespace?).

What decides this, I think, should be how we handle the code sharing.
My answer is to handle the code sharing with inheritance - it is a
natural, well-supported, flexible mechanism built into python for
doing exactly this sort of code sharing. All basis objects share some
features (zero(), __neq__()), and some share more (any graded basis
can sensibly act like a sequence, returning the ith basis function).
So it seems to me that inheritance is a sensible way to connect basis
objects and allow them to share code.

I've tried two different places to put the "meat" functions, directly
in polynomial classes, or in basis classes. The latter allows all
current polynomial representations to share a common type, which makes
implementing new representations a little tidier. Then the basis class
hierarchy allows code sharing between the "meat" functions too (e.g.
power, roots).

Could the basis objects share code in some other way? I don't really
see another good approach - for example I don't see how your
loose-functions-plus-class-maker would handle the issue of power(),
which is generic for some implementations but specialized in others. I
expect roots() to be the same.

>> ChebyshevBasis = generate_basis("ChebyshevBasis", mul=chebmul, ...)
>>
>
> Heh! Well it's does look like a compromise, or rather, design by committee
> ;) But it is something to think about.

You should have seen my first attempt - I tried to write a function
declare_function(globals(),"lag",LagrangeBasis) that would create
functional interfaces and add them to the module namespace, but I got
lost in the bowels of python's introspection system.

Joking aside, this really is applying your technique to my preferred
structure of polynomial object plus basis object. I think there's a
real need for a basis object, and if you want to be able to create one
simply by specifying a handful of low-level functions, this will do it
for you.

I should also say that, having implemented KroghInterpolator and
BarycentricInterpolator - polynomial objects in fact, albeit with a
different goal and very limited features - it's really not difficult
to write fast compiled code in cython that works in an object-oriented
way.

Anne