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

[Anne Archibald]

2009/10/13 Charles R Harris <charlesr.harris at gmail.com>:
>
>
> On Tue, Oct 13, 2009 at 10:12 AM, 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 8:14 AM, Pauli Virtanen <pav+sp at iki.fi> wrote:
>> [snip]
>> >> Anyway, you write the code, you call the shots :)
>> >>
>> >
>> > Heh. I am looking for feedback before taking the plunge.
>>
>> Wow, lots of ideas! But I think before we can evaluate them we need to
>> figure out what our requirements are. Here's a suggested list:
>>
>> (1) Simple. Simple to understand, simple to add new representations,
>> even for users outside scipy.
>> (2) Possible for users to subclass polynomials to provide polynomial
>> objects with more features.
>> (3) Polynomial objects should be able to detect and prevent operations
>> between polynomials in different bases when appropriate (for example,
>> you can substitute any kind of polynomial into a power-basis
>> polynomial, or a Lagrange basis polynomial into anything).
>> (4) Polynomial operations available as functions that operate on
>> coefficient arrays rather than polynomial objects.
>> (5) Possibly substantial data and methods attached to bases -
>> orthogonal polynomial roots and weights, lists of points for Lagrange
>> polynomials.
>>
>> Any others? Are all these really requirements?
>>
>
> I would add "massive code reuse".

I agree, which is why I used inheritance to connect bases and polynomial types.

>> Here are some key questions:
>>
>> * Do we want the information about what representation a polynomial is
>> in stored in the polynomial's class, in an explicit basis object the
>> polynomial contains, or both?
>
> I would say different representations, different classes, because different
> representations don't mix; you can't even add the coefficients together.
> This could be carried to an extreme, i.e., different scalings -> different
> classes, but I think that degree of specificity is too much. Also no class
> should need to know about any other class in order to work properly, it
> should stand independently.

The idea that they shouldn't mix is certainly correct, but having a
polynomial refer to its basis explicitly is another way to achieve
this; given the way python works, having separate classes doesn't even
allow you to avoid manual type checking.

class Foo:
    def __mul__(self, other):
        if not isinstance(other, Foo):
             raise ValueError

versus

class Bar:
    def __mul__(self, other):
        if not (isinstance(other, Bar) and other.basis==self.basis):
             raise ValueError

 or even, as with my current implementation:

class Bar:
    def __mul__(self, other):
        if not self.iscompatible(other):
             raise ValueError

I agree classes should not need to know about each other's internal
workings; this is part of requirement (2): we must be able to work
fine with polynomial representations added by users outside scipy.

>> * How should a user obtain a polynomial object?
>>
>
> I was thinking if you do 4), which is the fundamental requirement, then you
> should be able get an appropriate class on top of that without any extra
> work.

Well, my first attempt was a simple Polynomial object which was used
for polynomials in any representation; it delegated all the heavy
lifting to the basis object's (4)-style interface. Polynomial objects
were constructed in the usual way. Unfortunately, that means there's
no way for different bases to provide polynomial objects with more
methods. So I switched to using mybasis.polynomial(coefficients),
which lets the basis create a subclass if that's appropriate.

>> My current implementation fails only (4), and could be modified by
>> moving the code back into basis objects to satisfy (4). Then
>> polynomials of all current representations could be of class
>> Polynomial, with the basis information stored as an explicit basis
>> object. Users obtain polynomials as
>
> I really don't see any natural is-a relationship except for extending a
> class, say by adding weights. Which is to say I don't see any natural base
> class for all the different representations. I don't think inheritance is
> the right concept here.

Well, a GradedBasis is-a Basis, and a PowerBasis is-a GradedBasis, and
so on, but I think you were talking more about the polynomial objects
themselves. If polynomial objects delegate all their calculation to
the (4)-style interface, then is there any need for them to be
different types at all? After all, the delegation, compatibility, and
type-checking code is all common to all polynomial types.

In the Basis objects, there is room for massive code reuse - I have
generic companion matrix and root-finding code that works for any
polynomial representation that admits division, for example. Extending
coefficient arrays by padding with zeros works for any GradedBasis
(and in fact it's possible to write a division algorithm that works
for any GradedBasis). This takes natural advantage of a class
hierarchy among the bases.

As for polynomial objects, code reuse is most easily accomplished by
having them be exactly the same type; all the wrapper gubbins is the
same in all bases, apart from raising a few NotImplementedErrors.


Anne