[SciPy-dev] Generic polynomials class (was Re: Volunteer for Scipy Project)
Charles R Harris
charlesr.harris@gmail....
Thu Oct 15 10:55:38 CDT 2009
On Wed, Oct 14, 2009 at 10:53 PM, Anne Archibald
<peridot.faceted@gmail.com>wrote:
> 2009/10/14 Charles R Harris <charlesr.harris@gmail.com>:
> >
> > On Wed, Oct 14, 2009 at 7:09 PM, Anne Archibald <
> peridot.faceted@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).
>
> >> 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.
>
> 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. 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.)
>
>
> Anne
> _______________________________________________
> Scipy-dev mailing list
> Scipy-dev@scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-dev
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mail.scipy.org/pipermail/scipy-dev/attachments/20091015/4d71c169/attachment-0001.html
More information about the Scipy-dev
mailing list