Better adaptive plotting

Introduction

During Sage Days 9 I worked on three different projects. I finished two of them but I will only talk about one of them here.

The first one was to improve the rendering of plots in Sage. I implemented a better algorithm to insert additional points to the linear approximation of curves.

In all of the examples below, text looking like this:

some command
       
is meant as input to Sage, either in the notebook or at the command line and below it is the result it gives

Current algorithm

(until Sage 3.1.1)
The current plotting looks good for simple functions.
plot(sin(x), 0, 7) 
       
But if you try something wigglier, you get strange results.
plot(sin(1/x), -2, 2) 
       
plot(sin(x^4), -4, 4) 
       
You can increase plot_points to get better results (the default is 200).
plot(sin(1/x), -2, 2, plot_points=1000) 
       
plot(sin(1/x), -2, 2, plot_points=10000) 
       
plot(sin(x^4), -4, 4, plot_points=1000) 
       
plot(sin(x^4), -4, 4, plot_points=10000) 
       
And it looks pretty good after that, but the user shouldn't have to tweak parameters just to get a reasonable-looking plot.

All of this can be explained by the current algorithm that insert points in the rendering to adapt to the curve.
p = plot(x^2, -0.5, 1.5)
p += line([(0, 0), (1, 1)], rgbcolor='green')
p += point([(0, 0), (1, 1)], rgbcolor='red', pointsize=40)
p += point([(0.5, 0.25)], rgbcolor='purple', pointsize=40)
p += text('A', (-0.05, -0.1), rgbcolor='red')
p += text('B', (1.1, 1), rgbcolor='red')
p += text('C', (0.57, 0.2), rgbcolor='purple')
p.show(axes=False, aspect_ratio=1) 
       

If you have points A and B which are part of the sampled points on the curve, the approximation would be the line in green. The current algorithm would compute the difference between the y-coordinates of A and B. If it exceeds a certain threshold (which is 0.1 by default), then it adds the point C. It then repeats this procedure with C and B but not A and C. Otherwise it just proceeds to B and the next point. After having added some fixed number of points in this way (by default 10000), it just stops adding any more points.


Also if you start the plot at a singularity and you don't use a lot of points, you get rather strange results.
plot(sin(1/x), 0, 1, plot_points=5).show(ymin=0, ymax=1, xmin=0, xmax=1) 
       
Because currently, any point that cannot be evaluated is dropped from the plot.

The New Algorithm

(starting with Sage 3.1.2)
p = plot(x^2, -0.5, 1.5)
p += line([(0, 0), (1, 1)], rgbcolor='green')
p += point([(0, 0), (1, 1)], rgbcolor='red', pointsize=40)
p += point([(0.5, 0.5)], rgbcolor='red', pointsize=40)
p += point([(0.5, 0.25)], rgbcolor='purple', pointsize=40)
p += line([(0.5, 0.5), (0.5, 0.25)], rgbcolor='purple')
p += text('A', (-0.05, -0.1), rgbcolor='red')
p += text('B', (1.1, 1), rgbcolor='red')
p += text('C', (0.41, 0.54), rgbcolor='red')
p += text('D', (0.57, 0.2), rgbcolor='purple')
p.show(axes=False, aspect_ratio=1) 
       

If you have points A and B which are part of the sampled points on the curve, the approximation would be the line in green. The new algorithm would compute the point C and the y-coordinate difference between C and D. If it exceeds a certain threshold (which depends on the interval of x and the number of points), then it adds the point D. It then repeats recursivly this procedure with A, D and D, B. Otherwise it just proceeds to B and the next point. After having descended a certain number of recursion levels (5 by default) it stops recursing.


It still looks good for simple functions.
plot(sin, 0, 7) 
       
It's much better, with the default settings, for wiggly functions.
plot(sin(1/x), (-2, 2)) 
       
plot(sin(x^4), -4, 4) 
       
Also, as a special case for border points, if the value cannot be computed it tries other values that are closer and closer to the next point until it can compute the value.
plot(sin(1/x), 0, 1, plot_points=5) 
       

This means that if you start at a singularity, you still get an idea of the plot aspect. But since we used so little points this plot still looks bad. If we increase the adaptive_recursion factor, then it looks really good.

plot(sin(1/x), 0, 1, plot_points=5, adaptive_recursion=12)