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

[Charles R Harris]

On Tue, Oct 13, 2009 at 12:31 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 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<pav%2Bsp 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.
>
>
But those basis objects return a ton of NotImplementedErrors, which means
that somewhere down the line the methods have to be implemented. That was
one of the problems that I wanted to get around, because it leaves too much
code to be written. And makes having different implementations for the
methods, i.e., cython or Fortran, less transparent than just writing a set
of base functions (add, sub, etc.) in those languages. But I think we are
starting from somewhat different places. To me, the base functions were
primary and the classes conveniences built on top of them. That is more
flexible, IMHO, and doesn't tie one into a big system of class inheritance,
etc., which is one of the drawbacks I see in a lot in C++ projects. The GNU
scientific library would be an example of that where the pieces get tangled
up in the whole.

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.


> 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.
>
>
>
Chuck
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/scipy-dev/attachments/20091013/d2d9a425/attachment.html>