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

[Anne Archibald]

2009/10/13 Charles R Harris <charlesr.harris at gmail.com>:
>
>> I'm not sure I see why returning all those NotImplementedErrors is a
>> problem - yes, the code needs to be written to override the methods
>> that return NotImplementedError, but that code needs to be written no
>> matter how we set up an interface to it.
>
> It's not trivial to do a complete workup. My template class is 698 lines
> long, but it exists once in one place, and that does the trick for *all* the
> classes.

Aha. I think I see where the confusion is. In my proposed interface,
the coefficient-oriented functions would *be* the ones in the Basis
object. So, for example,

chebadd(c1,c2,interval) -> ChebyshevBasis(interval).add(c1, c2)

or, more likely, you'd have somewhere earlier

cheb = ChebyshevBasis(interval)

and then it's

chebadd(c1,c2,interval) -> cheb.add(c1,c2)

So the coefficient-oriented interface that we need to have is what
lives in the basis objects, and it is what the Polynomial object uses
to implement (e.g.) p1*p2. There would be no separate chebmul,
polymul, lagrangemul, ....

>> And having the base class
>> allows us to take any code that is the same between bases - like the
>> companion matrix code - and put it in the base class. So the only
>> "wasted" code here is the one-line wrappers that return
>> NotImplementedError, and these serve as documentation for users - any
>> such method is something they're going to have to implement for their
>> own polynomial representations. If you look at the three different
>> representations I have currently implemented, there isn't really much
>> repeated code.
>
> But it is also incomplete. And you have to have the separate functions
> anyway, so it is easier if you just need to write them and then you are
> done. If you want to change or fix the class interface, you don't have to
> touch the implementations of the functions. If you want to change the
> implementations, then you don't need to touch the interface. The two parts
> become mostly independent, which is generally a good thing.

The Basis object is effectively just a namespace holding the
collection of functions, and possibly any data they need (interval,
for example). So the separation you're describing is still present in
my proposed interface; it's the separation between Basis objects
(which contain the functions but no particular plumbing) and the
Polynomial objects (which are all plumbing).

>> > I think your method of getting the companion matrix is clever, and
>> > having
>> > definitions for zero, one, and x is a useful touch that should probably
>> > be
>> > part of any interface. I'm less convinced by the use of an empty list
>> > for
>> > zero.
>>
>> This choice is a natural one in some ways, because it means you never
>> normally have a leading coefficient of zero for polynomials of any
>> degree. It also means degree = len(coefficients)-1, apart from
>> polynomials that are represented in a degree-inflated manner. And the
>> places you need to special-case it are almost all places you would
>> need to special-case zero polynomials anyway. I realize it may not be
>> so convenient for a vector-of-coefficients interface, but I think even
>> there it's the natural approach.
>>
>
> But does it avoid special casing? Can you add it transparently to another
> array?
>
> In [1]: x = arange(3)
>
> In [2]: x + []
> ---------------------------------------------------------------------------
> ValueError                                Traceback (most recent call last)
>
> /home/charris/<ipython console> in <module>()
>
> ValueError: shape mismatch: objects cannot be broadcast to a single shape
>
> In [4]: x + [0]
> Out[4]: array([0, 1, 2])

This is a special case anyway - it's one of the very few instances in
which you can add coefficients the way you're describing. Normally you
have to extend one or the other so that they're the same length (and
in some bases this is a non-trivial operations, rather than simply
padding with zeros). If what you mean here is to add a scalar, there's
no reason to wrap it in a list;

In [3]: np.arange(3)+0
Out[3]: array([0, 1, 2])

In [4]: np.arange(3)+2
Out[4]: array([2, 3, 4])

If what you're adding is actually a vector of coefficients, you should
be surprised if it works:

In [5]: np.arange(3)+np.arange(2)
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)

/home/peridot/<ipython console> in <module>()

ValueError: shape mismatch: objects cannot be broadcast to a single shape

In [6]: np.arange(3)+np.arange(1)
Out[6]: array([0, 1, 2])

> I think it is a complication.

It's a question of trading one set of special cases for another. If
the zero polynomial is the only polynomial with zero leading
coefficient, that's going to force a fair number of special cases as
well. (Polynomials with zero leading coefficient that are generated as
a result of calculations are trouble of a different sort, since it's
entirely likely that the leading coefficient will be corrupted by
roundoff error and come out to something near zero, where the notion
of "near" is unclear.)

In any case, it's only really important for people who want to use the
coefficient-array interface; the polynomial class can hide this detail
under the rug. I'm not terribly attached to this representation.

Anne