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 < notifications@github.com> 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 https://github.com/tbreloff/Plots.jl/issues/621#issuecomment-268820984, 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.

[Armavica]

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

It works in Mathematica like a charm. 1661 points though.

discontinuities will always lead to max refinement

This is not necessarily a bad thing since it allows to bracket precisely the discontinuity point.

[KristofferC]

I would have thought the small intervals on the left side would have made the right side to start winning.

[KristofferC]

What does mathematica do for a square wave? How many points?

[KristofferC]

"subdividing a given interval at most MaxRecursion times". Maybe that is what is happening for Mathematica and sin(1/x)?

[Armavica]

67 points (discontinuity at 0.5).

[KristofferC]

It makes sense to limit the number of recursions because there is no point in having points too dense since they won't affect the figure. Should help with discontinuity problem as well. Will implement and report back.

[Armavica]

Interestingly Mathematica never evaluates the function at the boundaries of the interval, and treats them basically like discontinuities (it pushes towards them from time to time). This is indeed generally what you want…

[ChrisRackauckas]

Interestingly Mathematica never evaluates the function at the boundaries of the interval, and treats them basically like discontinuities (it pushes towards them from time to time). This is indeed generally what you want…

Is it always the same starting perturbation from the ends?

[Armavica]

It always starts at min+(max-min)/(n-1)*10^-6 (n number of initial points) for the minimum, and symmetrically for the max.

[KristofferC]

Version 4 https://gist.github.com/KristofferC/284559e8f856066ec3b922dadcd3df03, http://imgur.com/a/ZFTrf

Still some stuff left to do but getting better I think.

[KristofferC]

New (and probably final) version. Now, never evaluates the function exactly at the end points and uses some more sensible pruning of intervals. Before it over-refined regions it didn't need to just because the function was ill behaved in some other region.

https://gist.github.com/KristofferC/6b68b23527e6e9e668708764128c39cd results http://imgur.com/a/yV4HK

[Evizero]

very cool!

as a side note: also a nice function that was annoying to plot correctly is this one: f = x->abs(x)^.5, which is the LP distance loss with P < 1 (here concretely the blue one with P= 0.5), because of how steep it is near 0

[KristofferC]

Looks pretty good now?

[Armavica]

Other interesting functions to plot: lgamma on (-4, 8), tan on (-5,5) (for this one ylims needs to be specified manually, unless we implement some heuristics to detect the relevant y range). As a consequence of the nondeterministic plotting, the height of the peaks in lgamma varies from one plot to another, I think we should really consider if we want this to happen or not.

[KristofferC]

Yeah, manual ylims is also something worth considering. For the determinism, would it be good to just create a new RNG instance and give it a fixed seed. Then the noise is "random" but constant for the same input points.

[KristofferC]

Updated gist with using a custom rng object. This is lgamma: http://imgur.com/a/s4b4X

[ChrisRackauckas]

Wow, that looks really good.

[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 < notifications@github.com> wrote:

Wow, that looks really good.

— You are receiving this because you commented.

Reply to this email directly, view it on GitHub https://github.com/tbreloff/Plots.jl/issues/621#issuecomment-268996019, 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 notifications@github.com 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 https://github.com/tbreloff/Plots.jl/issues/621#issuecomment-268997354, or mute the thread https://github.com/notifications/unsubscribe-auth/AA492gZcaV4EmhUJQAbCzEIdJoE-GE-gks5rK9ntgaJpZM4LT8hR .