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

Wed Oct 14 19:38:09 EDT 2009

On Wed, Oct 14, 2009 at 12:43 PM, Anne Archibald
<peridot.faceted at gmail.com>wrote:

> 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.
>
>
I think we should look at Josef's suggestions of metaclasses also.

> >> 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.
>
>
Yes, I think the thing to keep in mind is what we want to do and then find
the best tool for the job.

> 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 don't like the idea of passing the basis into the instance constructor, I
think it would be preferable to pass it into the constructor of a class
creator. That is why Josef's suggestion of a metaclass rouses my curiosity.
Now if I can only figure out what metaclasses are...

> 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.
>
>
Yes. Although I expect the bottlenecks are few and those few can be pure C
(cython) code while the rest is just plain old python. Which cython will
compile also, of course ;)

Chuck
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/scipy-dev/attachments/20091014/536ae472/attachment.html>