[R] Polynomiographic function in R :-)

[François Pinard]

Hi, people.  Nothing too serious in this message.  Nevertheless, all
criticism or advice is welcome :-).

Yesterday, I went to a conference by Bahman Kalantari (Rutgers
University) about Polynomiography (the Fine Art and Science of
Visualizing Polynomials).  Since I'm starting my R learning, I decided
to try using it for computing some (any!) polynomiograph.  I was
surprised about how easy and quick it was to get some results.  Then I
thought it was interesting to extend the drawing to any function, and
not necessarily polynomials, yielding the function below:


polygraph <- function(expression, xrange=c(-1, 1), yrange=c(-1, 1),
                      points=200, steps=20, display=image)
{
    expression <- substitute(expression)
    variable <- all.vars(expression)
    stopifnot(length(variable) == 1)
    derivative <- D(expression, variable)
    name <- as.name(variable)
    expression <- substitute(name - expression / derivative)
    assign(variable, outer(seq(xrange[1], xrange[2], length=points),
                           seq(yrange[1], yrange[2], length=points) * 1i,
                           '+'))
    for (step in 1:steps) {
        display(Arg(eval(name)))
        assign(variable, eval(expression))
    }
}


which can be used this way, picking a function almost at random, say:


polygraph(x^3 - sqrt(x) - 1, points=300)


Here are a few random thoughts or remarks:

* Once fully converged, there should be only one colour per root.  Each
pixel colour shows towards which root would converge the chosen root
finding algorithm, starting at this particular point, or complex number.

* Another nice choice for `display' could be `filled.contour', yet it
computes more slowly.

* The successive plots (20 by default) show the progressive refinement
while finding equation roots, making a kind of animation.  One might
prefer moving the `display' call out of the loop, and show only the last
refinement.

* I did not know that root finding through Newton-Raphson could be
merely extended to complex numbers, fun to see that it works! :-)

* The conferencer told us that there a _lot_ of root finding algorithms,
and they may yield different styles of art.  I only picked the simplest
one to play with.  But you might do better!  (There are also many other
approaches than root finding for producing graphs out of polynomials.)

* Really, the one thing that most amused me in this experiment is how I
could use R for symbolically preparing the computation to do, without
resorting to parsing and deparsing (which I'm instinctively tempted to
avoid.)  I'm quite far from understanding all I should about functions,
expressions, calls and parse trees, but even knowing very little, it was
satisfying being able to rather quickly debug the above function.

* There are likely better ways than those I used.  For example, even if
unlikely, there might be clashes between the variables making up the
expression given, and local variables of the function.  I wonder if the
expression variable could have been more fully "abstracted".

* Vectorisation worked surpringly well on that problem, speed-wise.
However, because some regions of the plane converge faster than others
(use `display=plot' and such while calling `polygraph' to study this),
maybe they would be ways towards significant speed-ups.  But since it is
likely that one would loose a good part of "vectorisability" by doing
so, and add a lot of complexity (with unavoidable bugs in the process),
I wonder how worth it would be in practice.

* Given a matrix of complex results, they should ideally be turned
into N groups, each group being related to one of the N roots of the
equation.  I tried producing factors out of these results, but numerical
approximation made that non-practical.  I would guess that "clustering",
which I do not know, may be seen as a way to produce factors fuzzily.

* As a counter-measure to the above difficulty, I used `Arg()' as a way
to produce "levels" out of the results.  Could have used `Im()' instead.
It seems that `Mod()' and `Re()' are less productive. `image' is kind
enough to turn those levels into colours without any effort from me!


All in all, it is a fun way to explore R capabilities, and it also opens
up all kind of ideas to toy with! :-)

-- 
François Pinard   http://pinard.progiciels-bpi.ca