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

[Charles R Harris]

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

> 2009/10/14 Charles R Harris <charlesr.harris at gmail.com>:
> >
> > On Wed, Oct 14, 2009 at 7:09 PM, Anne Archibald <
> peridot.faceted at gmail.com>
> > wrote:
> >>
> >> I think passing it to the constructor is kind of a side issue (as I
> >> said, users are never supposed to call the constructor); if I
> >> understand correctly, what you're objecting to is the fact that every
> >> polynomial has a .basis attribute that points to the basis its
> >> coefficients represent.
> >>
> >
> > Wait a minute, let's get this straight. Users aren't supposed to use the
> > class object, but only instances generated from some primeval instance
> > created out of nothingness in the module? That sounds like a factory
> object,
> > but not of the traditional sort. Mathematicians may theoretically start
> with
> > the empty set and a couple of axioms, but most don't work that way ;) How
> > does repr work, since it is supposed to generate a string that can be run
> to
> > recreate the object? Doesn't it need to call the constructor? Same with
> > unpickling.
>
> Here's a typical use:
>
> b = ChebyshevBasis((0,1))
>
> X = b.polynomial([0,5,0.5])
>
> f = polyfit(basis=b, degree=1, x=[0,0.5,1], y=[1,2,1])
>
> Y = X**2 + f
>
> print Y([0,0.5,1])
>
> The constructor does indeed get called, but by the basis object
> (inside its polynomial method). It's still not possible to construct a
> polynomial without specifying its basis, of course; after all a
> coefficient array is meaningless without a basis in which to interpret
> it. So this may not resolve your objection to passing a basis to the
> constructor.
>
> >> Certainly every polynomial needs to know what basis it is in, and many
> >> users will also need to know. But it would be possible to create a
> >> class for every different basis and store the basis object in the
> >> class, so that a polynomial contains only its coefficients and a
> >> pointer to its class. Then you could extract a polynomial's basis by
> >> doing, say, mypoly.__class__.basis. Apart from the ugliness of doing
> >> this, I find it rather unpleasant to have potentially thousands of
> >> class objects floating around
> >
> > Well, identifying different classes by attaching a name attribute, in
> this
> > case a basis, is how some folks identified classes back in the
> paleolithic
> > period of C++, Borland comes to mind, IIRC, they used to attach an
> > identifier to plain old structures too. It works, sorta, but it is just a
> > lightweight version of having lots of different classes.  But it does
> seem
> > to me that your contruction is oriented towards Lagrange bases and
> > barycentric interpolation and the desire to have it also incorporate
> graded
> > polynomial basis and such has complicated it. It would be simpler to just
> > build a different classes for the different sorts of polynomials.
>
> Certainly my design is intended to accommodate the Lagrange
> representation. But it is also shaped by a desire to be able to deal
> sensibly with families of orthogonal polynomials, providing auxiliary
> information about them. There is also a certain amount of code that
> can be shared between the Lagrange basis polynomials and other
> polynomial implementations - all the type checking plumbing, for
> example. I'd be reluctant to pull them apart and duplicate that code.
> (I also rather like being able to use isinstance(x,Polynomial) but I
> realize that's not as pythonic as one might wish.)
>
> >> - remember PowerBasis(0) is different
> >> from PowerBasis(0.32). Also, you have to do some careful work to check
> >> polynomial compatibility: people can, and often will, create different
> >> basis objects that represent the same object. Either you need to
> >> implement a flyweight pattern for your basis objects (or classes -
> >> ugh), or you need to have a way of detecting when two polynomials with
> >> different classes are actually compatible (because the Basis objects
> >> used to manufacture the classes are equal but not the same object).
> >>
> > How do different classes end up with the same basis. And where do all the
> > different
> > basis come from if the user doesn't generate them and call the
> constructor?
>
> Well, in an ideal world, users would simply create one basis and use
> it for all the polynomials they produce in that basis. But it's
> already clear from writing the test code that we're going to have
> things like (in my syntax):
>
> p1 = PowerBasis(center=3).polynomial([1,1])
> ... many lines of code ...
> p2 = PowerBasis(center=3).polynomial([3,5,6])
>
> In this situation p1 and p2 have received different (i.e.
> id(A)!=id(B)) basis objects that are actually equal.
>
> As for different classes with the same basis, I think I misunderstood
> your proposal. Since your code requires bases to be specified
> completely by an implementation and an interval, and you store the
> interval in the polynomial object rather than the class, you never do
> end up with many classes. If I've understood you correctly, you do
> require the user to specify an interval any time they create a
> polynomial (inlcuding zero(), one(), and X()), and this is then
> checked for compatibility.
>
> I'm not sure I see how your code can be adapted to sensibly handle the
> power basis (specifying a center and maybe a scale) or some of the
> orthogonal polynomial bases (some of which need an endpoint and scale,
> or possibly also one or two parameters specifying the family of
> polynomials).
>
>
The interval defines an affine map from the specified interval to [-1,1],
and hence gives scale and center, it was just a guestion of choosing a
representation of that map. Since I almost always end up having to calculate
center and scale from two points, the interval approach seemed the most
natural. Thus plain old power series would be defined in the interval [-1,
1] by default, just like Chebyshev series. This breaks down a bit for the
Bernstein polynomials, which could just as easily be defined on [-1, 1] by
using the form (x + 1)^i * (x + 1)^(n - i), but that wouldn't be what folks
were expecting. This could be fixed by making the default interval a class
attribute.


> >> It just seems simpler to me to have the basis as a full-fledged
> >> attribute of the polynomial class. Given that polynomial classes are
> >> supposed to be immutable, that means it has to be passed in the
> >> constructor.
> >>
> >
> > I think is would be simpler to seek less generality. You have an abstract
> > Basis class, then a derived GradedBasis class, then a factory function to
> > generate a GeneratedBasis class derived from GradedBasis,  and that class
> > has to be instantiated somewhere. Then there is a Polynomial class that
> > accepts a basis (or basis class?). How is that simpler then a factory
> > function that just builds specific classes for specific instances of
> graded
> > polynomials out of a few defined functions and inherits only from object?
> > The setup you have looks a bit over designed.
>
> One reason it looks over-designed is that it only has three polynomial
> types; it should at least have one more for families of orthogonal
> polynomials, and perhaps another for Bernstein polynomials. (The
> orthogonal polynomials will inherit from GradedBasis but probably
> can't be produced using GeneratedBasis; the Bernstein basis is not
> graded.) The GeneratedBasis I would be inclined to do without, except
> that if people want to easily make polynomials given collections of
> coefficient-operating functions, it will let them do that. I should
> point out that, for example, there's no need to implement chebadd and
> chebsub, since those implementations are automatically provided by
> GradedBasis.
>
> Leaving aside the GeneratedBasis, I think it's fairly simple.
>
>
I think it would be interesting to see an implementation restricted to
Lagrange basis only and as simple as you could make it.


> From a user's point of view, they pick a polynomial representation
> they want to work with by creating an appropriate Basis object, then
> use that to create the polynomial, perhaps directly from coefficients
> (with basis.polynomial()) or roots (basis.from_roots()), perhaps using
> existing ones (basis.X(), basis.convert(other_polynomial)), or by
> using some polynomial-returning function (e.g. poly.polyfit(basis,
> ...)). They can then use that polynomial as an immutable type.
>
> From an implementer's point of view, adding a new polynomial type is
> simple: they define a new basis class inheriting from the most similar
> existing basis class, and implement a few operations on coefficient
> arrays. Addition and power they can leave out; if they don't care
> about division they can leave it out; if they implement division or
> companion matrices they get root-finding for free. They'll probably
> want to write a routine to convert arbitrary polynomials to their
> representation; the other direction is already implemented.
>
> If they've got low-level functions like chebmul to work with, they can
> simply write MyBasis = GeneratedBasis(...functions...) and lo and
> behold they've got a reasonably functional polynomial class. Or if we
> ditch that as too messy, they write a class with some one- and
> two-line wrappers.
>
> If you want to compare our two designs, mine uses a single class for
> all polynomials and uses the different classes of Basis to provide
> different implementations, while yours puts the code directly in the
> polynomial classes. Mine uses inheritance to share code between
> representations, while yours shares code by writing it into the class
> generator function.
>
> I think that to the extent mine is more complicated, it's because it
> wants to allow more different implementations of polynomial classes.
> Yours is very specialized to orthogonal polynomials defined on a
> finite interval. You'll need to generalize it somewhat to handle
> power-basis polynomials at least. I find the class-generating code
> fairly opaque, and I'm not sure how well it will interact with
> debugging and introspection. If these aren't a problem, and we really
> don't need more polynomial types, then maybe it does make sense to go
> with yours as the simpler.
>
> I'd argue, though, that a general, efficient, numerically-stable
> polynomial class could replace KroghInterpolator and
> BarycentricInterpolator and serve to make PiecewisePolynomial more
> powerful, for example.


I agree completely on this point. That is one reason I think a "simplest"
class for Lagrange bases would be interesting to look at.


> The orthogonal polynomial code in scipy.special
> could return polynomial objects, ideally in their own basis, that
> allowed fast evaluation, integration, differentiation, and whose basis
> objects carried information about the roots and weights needed for
> Gaussian quadrature. (In fact, those basis objects themselves could
> live in scipy.special and serve to produce the orthogonal polynomials
> as needed.)
>
>
All you need to do for the Gaussian quadrature is to carry along a weight
function.  But the quadrature module already has that bit of infrastructure,
it uses the standard companion matrix defined by the recursion along with
the weight function. My feeling was that having that as part of the basic
polynomial class was starting to overburden the class with accessories.

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