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

[Charles R Harris]

On Fri, Oct 16, 2009 at 1:58 PM, Anne Archibald
<peridot.faceted at gmail.com>wrote:

> 2009/10/15 Fernando Perez <fperez.net at gmail.com>:
> > On Thu, Oct 15, 2009 at 6:36 PM, Charles R Harris
> > <charlesr.harris at gmail.com> wrote:
> >>
> >> How did you build the analysis objects? Part of what we have been
> discussing
> >> it the easiest way to expose a consistent interface while reusing the
> >> maximum amount of code.
> >
> > See
> >
> > http://github.com/fperez/nitime/blob/master/nitime/timeseries.py
> >
> > starting around line 504.  The analyzers store a reference to a time
> > series object and then expose properties (that are replaced in-place
> > by a static attribute once computed, as explained in the paper) that
> > make the library calls.  All the actual numerics are in
> >
> > http://github.com/fperez/nitime/blob/master/nitime/algorithms.py
> >
> > I'm not saying this will exactly solve your current design issues, but
> > it may provide some useful ideas to think about.
>
> Well, I think I'm convinced of the necessity for a low-level
> procedural interface, so I've set up a draft one in git. There's some
> question of just how bare-bones it should really be - for example
> should it include an add function, since that's just using the
> coefficient-extension function then adding? What about multiplying,
> which is the same thing, though this is not obvious?
>
>
It should definitely include add. For my personal use the following tend to
be most important:

Stuff I use:

   - evaluate
   - differentiate
   - integrate
   - leastsquares
   - zeros

Stuff I don't use much but should be there:


   - add
   - subtract
   - multiply
   - divide

Stuff for completeness:

   - powers

Conveniences?

   - weight function
   - one
   - zero
   - X




> As it stands I've only included implementations of functions with
> nontrivial algorithms.


Include them all, even the trivial ones. It may look silly, but more often
than not it will save hassles down the line. Remember, you aren't doing this
for you and *your* classes, you are doing it for someone else and *their*
classes. Why should they have to worry about how you do addition? If nothing
else the meaning of the inputs will be documented.


> This has the advantage that, realistically,
> these functions will not be shared between different polynomial
> representations. If I start including more straight-forward functions,
> I will need to figure out how to share the implementation of "add"
> between different representations if I want to avoid it being
> boilerplate code for all representations.
>
>
Write a polyutils module and define:

my_utterly_trivial_addition_function_that_bores_me =
common_utterly_trivial_addition_function

If we're implementing the heavy lifting in a procedural interface,
> then the object implementations will just be plumbing (as I see your
> analyzer objects are). So on the one hand there's not much code to be
> shared, but on the other it's all boilerplate that would benefit from
> code sharing.
>
>
This is where the polynomial code differs from the nipy code: they have
analyser objects with different interfaces that use the same basic functions
on different datasets. We want to have a common "analyser" interface that
uses different basic function.


> I think, though, that the next important decision has nothing direct
> to do with code sharing; we need to settle on what polynomial objects
> will look like to users.
>
> Some things are clear:
>
> * Polynomial objects will be immutable.
> * Polynomial objects will support + - * / % ** divmod and __call__
> * Polynomial objects will be able to report their roots
> * Polynomial objects will be able to return polynomials representing
> their derivatives and antiderivatives
> * Some polynomial objects will have extra features, like the ability
> to trim small coefficients from Chebyshev polynomials
> * We need some method of converting from one representation to another
>

Since the Chebyshev class can operate with objects, evaluating the Chebyshev
series for the "X" polynomial object does the conversion.  This works for
all the graded polynomial classes and might could be made to work with
Lagrange basis. Hopefully, that is something that isn't used a whole lot
anyway. Efficient evaluation for numeric arguments need to be reasonably
efficient, evaluation for objects less so.


> * We should have convenience objects zero, one, and X, to build
> polynomials.
> * There should be a least-squares polynomial fitting function.
>

+1

Maybe a standard interval also?


> * Some polynomials will have functions that generate them (e.g.
> Chebyshev series for a function)
>

Isn't this a special case of fitting?


> * We need at least power basis, Chebyshev, and Lagrange polynomials.
> * Ultimately the polynomial class should be able to replace
> KroghInterpolator and BarycentricInterpolator and be stably returned
> from the orthogonal polynomial routines
>
>
We probably want trigonometric polynomial class(es) also. I use barycentric
forms of these for the complex remez algorithm. They are just streamlined
versions of a lagrange basis with the sample points on the unit circle in
the complex plane or the two-sheeted Riemann surface equivalent.


> Less clear:
>
> * Should polynomials store a degree attribute or should it be inferred
> from the size of their coefficient array? (Sometimes it's handy to
> work with degree-inflated polynomials, but deflating their degree is a
> nontrivial operation.)
>

I choose the latter because it is more dynamic. The coefficients are public
attributes and if some fool wants to assign a different array to them I
figure thing should still work. Computing it is unlikely to be a performance
bottleneck in any case.

* Should we try to accommodate complex variables? Complex coefficients?
>

Complex coefficients, definitely. I might not use them, but no doubt someone
will want them. If there are stability problems, that is part of the game
and maybe the user can offer improvements. Caveat emptor.


> * Should we try to accommodate vector values?
>

If you mean evaluate at a vector of points, then yes, it is the only way to
overcome function calling overhead if you need a ton of values.


> * Should polynomials have a method to evaluate derivatives and
> antiderivatives directly without constructing the
> derivative/antiderivative object first?
>

I think that is unnecessary.


> * Do we need Lagrange-Hermite (some derivatives specified as well as
> values) polynomials? (KroghInterpolator implements these)
>

Don't know


> * How should code that works for any polynomial type (e.g.
> multiplication by a Lagrange polynomial) identify that an object it
> has been given is some sort of polynomial? Callable and presence of a
> degree attribute? Subclass of Polynomial?
>

I don't think the classes need to interoperate, in fact I would discourage
it because it complicates everything for a need I don't see. But if there
are such functions, then deriving everything from a dummy Polynomial class
would be the best approach IMHO. We could do that anyway to leave the door
open to future changes.

class Polynomial: pass



> * How should users implement their own polynomial types so they
> interoperate with scipy's built-in ones?
>
>
What do you mean by interoperate?


> Still debated:
>
> * I think we need a first-class object to represent the basis used to
> express a polynomial. Charles Harris disagrees; he wants to encode
> that information in the class of the polynomial plus an interval
> stored in the polynomial. Some different scheme would necessarily be
> used for Lagrange polynomials.
>

Things are encoded in 1) an array of coefficients 2) an interval, and 3) the
base functions and 4) the class name. In the Nipy model, 1) and 2) would be
the data containers.


> * How should the zero polynomial be represented? I think its
> coefficient array should be np.zeros(0); Charles Harris prefers
> np.zeros(1).
>
>
You know where I stand ;) For a Lagrange basis of indeterminate degree I
think is should be an array of zeros of the appropriate dimension.

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