Plot Recipe for Automatic Adaptive Plotting of Functions

[ChrisRackauckas]

It would be nice to be able to plot(f) for f(t) some function which returns scalars (or more). Mathematica handles this with a sophisticated algorithm. It's only mentioned on this page:

https://reference.wolfram.com/language/ref/Plot.html

and described as

Plot initially evaluates f at a number of equally spaced sample points specified by PlotPoints. Then it uses an adaptive algorithm to choose additional sample points, subdividing a given interval at most MaxRecursion times.

It seems like the algorithm does the following (or at least something like it):

1. Take an initial number of points (lets say 100)
2. Take two points skipping over one, do a linear interpolation
3. Check the difference between the linear interpolation and the actual value (the skipped point). If it's large enough, subdivide the interval.
4. Recurse 3 on the points you are getting until max_recursion or it hits the tolerance.
5. Repeat 2 until you're point of point triples.

This form of adaptive plotting would make it easy to plot a function without knowing the details. This is nice for when you don't actually know what the function looks like until plotting it!

@Armavica

[KristofferC]

I was bored so I did the simplest thing I could: http://imgur.com/a/SksX3

Script is here: https://gist.github.com/KristofferC/b13f2375e4b7596db7b2a82237151709

[KristofferC]

I'm refining 50% of the intervals each iteration because the total cost of the function evaluation should then be just an extra factor of two compared to the last evaluation. (2^0 + 2^1 + ... 2^n) = 2 * 2^n

[tbreloff]

Awesome. We could make this be the new recipe for single-function calls

…

On Thu, Dec 22, 2016 at 8:20 AM Kristoffer Carlsson < ***@***.***> wrote: I'm refining 50% of the intervals each iteration because the total cost of the function evaluation should then be just a factor of two compared to the last evaluation. (2^0 + 2^1 + ... 2^n) = 2 * 2^n — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <#621 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AA492krgGldAbPws9rUbx07UfaY1W0S3ks5rKpUegaJpZM4LT8hR> .

[KristofferC]

But you still need to give start and end points, I guess?

Also, how would you differ between plot(f, [0.0, 0.1]) where [0.0, 0.1] could either be the interval to plot on, or two individual points?

[KristofferC]

This implementation is of course also susceptible to aliasing, for example if you create a sine wave that is zero at two adjacent initial "seed" points. Maybe sprinkling some random points in each iteration could help. There must be tons of literature for this though. Not that I am personally not so interested in reading it :P

[Armavica]

Without looking at your posts I wrote exactly the same algorithm (also refining half the points). Mathematica seems to deal with the aliasing issue with a nontrivial periodic (?) pattern in the initial sampling: https://imgur.com/a/UC1vv (I plotted the differences between consecutive initial sampling points divided by the average interval length, so for a regular sampling everything would be 1).

[ChrisRackauckas]

Oh wow, you guys really ran with this!

This implementation is of course also susceptible to aliasing, for example if you create a sine wave that is zero at two adjacent initial "seed" points. Maybe sprinkling some random points in each iteration could help. There must be tons of literature for this though. Not that I am personally not so interested in reading it :P

I presume there's some randomness at each step which gives new points. If you look at the picture from Mathematica at the top, you see that the points don't where the midpoints would be, they seem slightly off. I was wondering why this would be the case, and aliasing would be a really good explanation. I think throwing in a little bit of noise in the initial selection of points, and a little bit of noise in each interval subdivision (amount of noise as a percentage of the interval size?), and it should be good should be good enough.

Also, how would you differ between plot(f, [0.0, 0.1]) where [0.0, 0.1] could either be the interval to plot on, or two individual points?

What about plot(f, (0.0, 0.1)) as you had it? As a tuple, that means you can also dispatch on the fact that it has to be only two points: start and end. Or you can extend this a little by, if given more points, force a plot at those points. Like plot(f, (0.0,0.5, 0.1)) will have an initial point at 0.5. Though probably the parameters used should just be exposed as keyword args.

[Armavica]

If you look at the picture from Mathematica at the top, you see that the points don't where the midpoints would be, they seem slightly off.

Actually running Mathematica code suggests that the refinements are always done in the middle of two consecutive points. I suspect that this picture is just an "artist view" and not to be taken too seriously.

That said, we are not a Mathematica clone and I don't see any problem with some randomness here (except that this would make plotting nondeterministic, for "nice" functions it shouldn't be too bad but maybe we don't want that).

[Armavica]

I am also thinking that we could easily extend this to parametric plots.

[ChrisRackauckas]

I think, going along with the idea I had about giving more timepoints, instead it would be really nice to be able to declare known points of discontinuity, and ensure that the initial points straddle it dis+epsilon/2 and dis-epsilon/2 away. Being able to declare discontinuities might reduce the number of points in some cases? But from @KristofferC 's pictures I don't see this as a real problem.

[Armavica]

Unless I am mistaken each discontinuity will generate at most 2*maxrecursions extra points, right? (4*maxrecursions in @KristofferC's version). It doesn't seem to be a huge cost.

[KristofferC]

I will fix something better than last one. Just gotta dress up the Christmas tree zzz.

[KristofferC]

Ok, new version at: https://gist.github.com/KristofferC/ce5169000b29cc6b15588dc4e133ee5c with results: http://imgur.com/a/iBAi2.

This one should just estimate the second derivative in each point and also wiggle the points a bit to try to remove degenerate cases. Not sure how to decide which one is best though...

Also, like Tom said, this could be used for simply plot(f) where the limits are taken from the xlim of the plot. If no limits are specified then the default ones are simply used.

[KristofferC]

I think it is important not to have too fine grid in the start as well. 50 points is probably too much.

[ChrisRackauckas]

It seems like this second one is safer, so I'd assume it is 👍 . The only way to really know is to let it in the wild for a year and see what crazy things end up causing a problem.

50 points is probably too much.

And could that be a kwarg? Along with rand_factor? I think some way of using an error tolerance to stop early is really crucial though, because your plots do have a lot of points in them.

[KristofferC]

The function right now has a max_points kwarg but not in inital_n_points kwarg. I don't think the rand_factor need to be user definable. Do you have a suggestion for an early exit criteria?

[KristofferC]

I think multiplying the error with the interval should be done as well. Basically integrating the second derivative over the interval. I updated the gist.

[ChrisRackauckas]

Exit early if sum(abs2,err)<tol?

[ChrisRackauckas]

I think you can just generate err once before the loop, and resize! in the loop. Maybe that could reduce some allocations due to how Julia caches it? I'm not sure how much that would really help though.

(same for fs)

[KristofferC]

Yeah it is not optimized but to be honest, it is not really needed. It goes fast enough for simple functions and for complicated functions, time would be spent in evaluating the function.

I don't think sum(abs2, E) is so good because it will be unit dependent and also dependent on the number of points.

[ChrisRackauckas]

Then what about something like abstol + reltol of ODE solvers?

est = sum(abs2,err./(abstol + max(f.(xs))*reltol))/length(err)

Then stop when est<1? abstol would need units to be the same as err. Maybe it's easiest to have only a relative tolerance here. For plotting, you can really only "see" relative error anyways.

Edit: that max isn't the error, it's the values

[KristofferC]

My use of err in the code is not good terminology. An error should go to zero as the number of points goes to infinity. It is not actually an error, it is an estimation of the integral of the second derivative over a small interval and I use that to chose what intervals to refine.

In adaptive finite element, one way of estimating the error is to see what happens to the solution after a hierarchial p- or h-refinement where you project the solution of the coarse mesh onto the fine mesh (or higher polynomial mesh). I will try something like that.

[ChrisRackauckas]

I see. I thought it was a comparison. Then yes, a dual mesh comparison for an error is the way to go. Though you could estimate errors without resorting to a second mesh by checking the inbetween points against a linear interpolation.

[KristofferC]

Updated with error estimator: https://gist.github.com/KristofferC/4b999e8b376d5073f7d4d8f69cf3352e.

Results at: http://imgur.com/a/15WsA

For linear function we break straight away. Other function breaks at "decent" point density I think. Feel free to test with your own functions.

[Evizero]

Fancy!

[Armavica]

I think that the real test is sin(1/x). Currently it does show some trouble to plot it between 0 and 0.1 (or rather 1e-6 and 0.1 since it does not like 1/0). Even pushing the max_point limit to 1000 a lot of points get wasted on the left part. Maybe we could detect that a region is overcrowded and that it is useless to invest any more points in it.

[KristofferC]

Yeah, that is the evil one :P Does other software handle it "well"?.

[Evizero]

I am getting an error when trying a zero one loss function on 0.5, but that could very well be my fault somehow in my blind ignorance of just calling it with no careful thought whatsoever.

Example plot of the function:

julia> f = x -> x > 0 ? 0 : 1
(::#18) (generic function with 1 method)

julia> mm = (-2, 2)
(-2,2)

julia> xs = refine_grid(f, mm)
ERROR: InexactError()
 in setindex!(::Array{Int64,1}, ::Float64, ::Int64) at ./array.jl:415
 in #refine_grid#17(::Int64, ::Float64, ::Function, ::Function, ::Tuple{Int64,Int64}) at ./REPL[87]:67
 in refine_grid(::Function, ::Tuple{Int64,Int64}) at ./REPL[87]:3

[KristofferC]

It was just an array that got turned into Int from your function and then the interpolation failed. I updated script.

[KristofferC]

Right now I remove the top 5% of the maximum errors for the termination criteria in order for the error not to blow up at discontinuities. No idea if it is a good idea but otherwise, discontinuities will always lead to max refinement.

[tbreloff]

Great stuff! Do you know where to add this code in Plots? (Look for the single function and function with endpoints recipes, I think in recipes.jl or series.jl?)

…

On Fri, Dec 23, 2016 at 7:15 AM Christopher Rackauckas < ***@***.***> wrote: Wow, that looks really good. — You are receiving this because you commented. Reply to this email directly, view it on GitHub <#621 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AA492gpi-6OFq_EHsg3h6PdYT3MjSlElks5rK9eMgaJpZM4LT8hR> .

[Evizero]

Maybe this would fit better into PlotUtils? so that freeloaders such as the UnicodePlots guy (i.e. me) can use it as well (if he ever gets around to his refactor that is). Anyhow, big fan either way

[tbreloff]

I'd be cool with adding the core logic to PlotUtils and just calling it from the Plots recipe

…

On Fri, Dec 23, 2016 at 7:25 AM Christof Stocker ***@***.***> wrote: Maybe this would fit better into PlotUtils? so that freeloaders such as the UnicodePlots guy, can use it as well (if he ever gets around to his refactor that is). Anyhow, big fan either way — You are receiving this because you commented. Reply to this email directly, view it on GitHub <#621 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AA492gZcaV4EmhUJQAbCzEIdJoE-GE-gks5rK9ntgaJpZM4LT8hR> .

[KristofferC]

Ref JuliaPlots/PlotUtils.jl#6.

Regarding explicit y-limits I guess it could be solved by returning many sub intervals where the function is inside the y-limits on all the intervals. However, I don't think I can allow myself to work much more on this (not that it isn't fun).

[KristofferC]

Another way of dealing with y-limit is simply not to refine intervals where the function is outside the y-limits.

[tbreloff]

This is great work... thanks!

using Plots; scatter(tan)

[soegaard]

FWIW I have collected a small set of horror functions to test plotting algorithms:

https://github.com/soegaard/bracket/blob/master/plotting/adaptive-plotting.rkt#L225

A lot of them are from
www.mines-paristech.fr/~kajler/pub/9507-RAOBNK-issac95.ps

[mkborregaard]

Have you tried plotting them with Plots?

[KristofferC]

FWIW, I am sure the stuff I implemented is not at all state of the art, quality wise. It was just me applying what sort of makes sense. It would be interesting to see how it works on these horror plots.

[soegaard]

No, I haven't tried Plots - I found this discussion through a comment on https://scicomp.stackexchange.com/q/2377/1541

[mkborregaard]

Are you familiar with Julia? Could you describe the racket syntax clearly enough for someone here to port them to Julia? It could be fun to see.

[soegaard]

I'm not familar with Julia. The Racket syntax is just prefix notation. Instead of a+b we write (+ a b).
So

     (+  (tan (+ (expt x 3) (- x) 1)) 
             (/ (+ x (* 3 (exp x)))))) 
    (interval xmin=-2 to xmax=2)

is tan( x^3 -x + 1) + 1 / (+ x 3 e^x)

Here (expt a b) means a^b and (exp a) means e^a. Finally the one argument version of / is the same as 1/x.

Other functions on the list:

sin(1/x)
sin(x^4)
sin(300x)
1 + x^2 + 0.0125 log( abs(1- 3*(x-1)))
sin(e^x)
1/sin(x)
sin(x)/x
log(1+sin(cos(x)))

These are simpler. They check specific properties:

  1/x            ; check division by zero and singularity at x=0
  tan(x)       ; check multiple singularities
   -1/x         ; check symmetry
  1/abs(x)   ; check connected components of same sign
  1E6 * x     ; check that lines with high slopes are drawn
   1E50 * x ; check that lines with high slopes are drawn

-2. 2. -15. 15.)

[mkborregaard]

Here's the first function on the list

f(x) = sin(1/x)
using Plots
plot(f)

plot(f, -2, 2)

plot(f, linspace(-2,2,1000))

[KristofferC]

That looks awful lol.

[KristofferC]

I thought I tested sin(1/x) specifically...

[pkofod]

I'm not familar with Julia. The Racket syntax is just prefix notation. Instead of a+b we write (+ a b).
So

shouldn't be too hard too parse that as Julia, just need to add commas and put the operator outside the parentheses

julia> +(2, /(4, 3))
3.333333333333333

[KristofferC]

I think there is a difference with how max deals with NaNs in 0.6 or something...

[KristofferC]

JuliaPlots/PlotUtils.jl#20

[mkborregaard]

@pkofod great that it's that simple! So do all the functions plot as prettily as sin(1/x) after @KristofferC 's update?

[KristofferC]

Update what?

[KristofferC]

Oh lol nvm.

[pkofod]

@pkofod great that it's that simple! So do all the functions plot as prettily as sin(1/x) after @KristofferC 's update?

I don't know... I've banned myself from using Pkg.anything until I've presented today!

[KristofferC]

What is really needed is to take the ylimits into account when doing the partitioning.