[WIP] Add a very basic Makie extension

[asinghvi17]

Checklist

• Appropriate tests were added
• Any code changes were done in a way that does not break public API
• All documentation related to code changes were updated
• The new code follows the
  contributor guidelines, in particular the SciML Style Guide and
  COLPRAC.
• Any new documentation only uses public API

Additional context

I thought I'd update the code in #427 for the latest Makie, and while I was at it made it an extension as well.

This doesn't work fully yet, in that:

1. It is only hooked up to AbstractTimeseriesSolutions.
2. It is currently only able to plot line plots.
3. It currently does not accept theme keywords passed to plot, like linewidth. Setting that theme externally still works, though.

[asinghvi17]

Ok, made some decent progress here! The recipes are now to the point that they can plot any DiffEq solution to any point-based plot type, depending on what the user provides. Color cycling also comes naturally and is affected by theme.

On the Makie level, it's currently not supported to pass themes from the plot invocation to subsidiary recipes unless they are specified and consumed, so you can't adjust things like linewidth in the plot call (though you can do that in the global theme, which does trickle down). Once that's solved on the Makie end it should be a relatively simple update here.

I also added basic support for ensembles, which seem to be the only things which had special treatment in SciMLBase's Plots recipes so far.

Here are some examples

Code here

using DifferentialEquations, Makie
using CairoMakie

function lorenz(du,u,p,t)
    du[1] = p[1]*(u[2]-u[1])
    du[2] = u[1]*(p[2]-u[3]) - u[2]
    du[3] = u[1]*u[2] - p[3]*u[3]
end

u0 = [1., 5., 10.]
tspan = (0., 100.)
p = (10.0,28.0,8/3)
prob = ODEProblem(lorenz, u0, tspan,p)
sol = solve(prob)

f, a, p = plot(sol)
axislegend(a)
f

f = Figure()
topgl = GridLayout(f[1, 1])
botgl = GridLayout(f[2, 1])
xyzt, _p = plot(topgl[1, 1:2], sol, plotdensity=10000,linewidth=1.5)
xy, _p = plot(topgl[1, 3], sol, plotdensity=10000, idxs=(1,2))
xz, _p = plot(botgl[1, 1], sol, plotdensity=10000, idxs=(1,3))
yz, _p = plot(botgl[1, 2], sol, plotdensity=10000, idxs=(2,3))
xyz, _p = plot(botgl[1, 3], sol, plotdensity=10000, idxs=(1,2,3))

axislegend.((xyzt, xy, xz, yz, xyz))
f

_f(x,y,z) = (sqrt(x^2+y^2+z^2),x)
plot(sol,idxs=(_f,1,2,3))
sol

f = Figure()
lines(f[1, 1], sol)
scatter(f[2, 1], sol; denseplot = false)
scatterlines(f[3, 1], sol; denseplot = false, plotdensity = 1)
f

function fx(du, u, p, t)
    du[1] = p[1] * u[1] - p[2] * u[1] * u[2]
    du[2] = -3 * u[2] + u[1] * u[2]
end

function gx(du, u, p, t)
    du[1] = p[3] * u[1]
    du[2] = p[4] * u[2]
end

p = [1.5, 1.0, 0.1, 0.1]
prob = SDEProblem(fx, gx, [1.0, 1.0], (0.0, 10.0), p)

prob_func = let p = p
    (prob, i, repeat) -> begin
        x = 0.3rand(2)
        remake(prob, p = [p[1], p[2], x[1], x[2]])
    end
end

ensemble_prob = EnsembleProblem(prob, prob_func = prob_func)
sim = solve(ensemble_prob, SRIW1(), trajectories = 10)

plot(sim)
plot(sim; idxs = (0, 1))

summ = EnsembleSummary(sim, 0:0.1:10)

plot(summ)

A basic plot

Support for a random recipe which is defined as consuming point-based input

Theming :)

Ensemble

Ensemble with errorbars

[ChrisRackauckas]

I also added basic support for ensembles, which seem to be the only things which had special treatment in SciMLBase's Plots recipes so far.

Integrators have a recipe as well

[ChrisRackauckas]

I think this looks good?

[codecov]

Codecov Report

Attention: 134 lines in your changes are missing coverage. Please review.

Comparison is base (48771ae) 40.16% compared to head (0636e40) 38.91%.
Report is 2 commits behind head on master.

❗ Current head 0636e40 differs from pull request most recent head 57890b3. Consider uploading reports for the commit 57890b3 to get more accurate results

Files	Patch %	Lines
ext/SciMLBaseMakieExt.jl	0.00%	133 Missing ⚠️
src/ensemble/ensemble_solutions.jl	0.00%	1 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #611      +/-   ##
==========================================
- Coverage   40.16%   38.91%   -1.26%     
==========================================
  Files          54       55       +1     
  Lines        4133     4266     +133     
==========================================
  Hits         1660     1660              
- Misses       2473     2606     +133     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[asinghvi17]

It should work for most simple usecases - more advanced ones will probably need to be deconstructed manually, but that's not particularly hard...

[ChrisRackauckas]

That's fine, it's a start. We can get it in there and let it keep improving.

[henry2004y]

I have a question about the recipe, but I am not sure if this is the right place to ask @asinghvi17 . Say if I want to scale the axis (e.g. time) by a factor of 10 for plotting, how should I do that? Makie provides a function scale!; it works for a simple case like

f = Figure()
ax = Axis(f[1, 1], aspect=1)
l0 = lines!(ax, [0, 1], [0, 1])
scale!(l0, 10, 1)

but not for the DiffEq recipe:

using OrdinaryDiffEq
using GLMakie

function lorenz(du,u,p,t)
    du[1] = p[1]*(u[2]-u[1])
    du[2] = u[1]*(p[2]-u[3]) - u[2]
    du[3] = u[1]*u[2] - p[3]*u[3]
end

u0 = [1., 5., 10.]
tspan = (0., 100.)
p = (10.0,28.0,8/3)
prob = ODEProblem(lorenz, u0, tspan,p)
sol = solve(prob, Tsit5())

f, a, p = plot(sol)
axislegend(a)
# This is not working!
scale!(p, 10, 1)

f

[asinghvi17]

Ah! I think this is another issue uncovered - the API I used to do this is still pretty new!

You can of course scale the plots individually, by inspecting e.g. ax.scene.plots, but this seems like an issue which has to be fixed in Makie. The primary thing that has to be solved is that a plot's Transformation should be passed to its children in the event that a plotlist is being generated.

Would you mind filing an issue on the Makie repo, @henry2004y?

[henry2004y]

Sure! Hopefully there will be a quick fix. Thanks for the new recipe!

[henry2004y]

And another maybe related issue is that the new recipe does not work with plotting attributes (e.g. linestyle, linewidth)? Is this an issue for Makie or the recipe @asinghvi17 ?

MWE

using OrdinaryDiffEq
using GLMakie

function lorenz(du,u,p,t)
    du[1] = p[1]*(u[2]-u[1])
    du[2] = u[1]*(p[2]-u[3]) - u[2]
    du[3] = u[1]*u[2] - p[3]*u[3]
end

u0 = [1., 5., 10.]
tspan = (0., 100.)
p = (10.0,28.0,8/3)
prob = ODEProblem(lorenz, u0, tspan,p)
sol = solve(prob, Tsit5())
f, a, p = lines(sol, linestyle=:dot)

[asinghvi17]

Same issue unfortunately :(

[henry2004y]

Same issue unfortunately :(

Is there a temporary workaround for attributes before it's properly fixed in Makie?

[henry2004y]

Honestly I am now pretty frustrated about the new SpecAPI introduced in Makie v0.20. The old approaches for creating recipes did not limit it capability in so many ways.

Here is yet another issue with PlotSpec: it does not work (yet) with lift:

MWE

using OrdinaryDiffEq
using GLMakie

function lorenz(du,u,p,t)
    du[1] = p[1]*(u[2]-u[1])
    du[2] = u[1]*(p[2]-u[3]) - u[2]
    du[3] = u[1]*u[2] - p[3]*u[3]
end

u0 = [1., 5., 10.]
tspan = (0., 100.)
p = (10.0,28.0,8/3)
prob = ODEProblem(lorenz, u0, tspan,p)
sol = solve(prob, Tsit5())

fig = Figure()

ax = Axis(fig[1, 1])

t1, t2 = 0.0, 100.0
step = (t2 - t1) / 1e3
time_format = t2 < 2e-3 ? "{:.3e} s" : "{:.3f} s"
# the minimum of the range cannot be t1.
sg = SliderGrid(fig[2, 1],
    (label="Time",
     range=range(t1+step, stop=t2, step=step),
     format=time_format,
     startvalue=t2))
# convert tspan to an Observable
tspan = lift(sg.sliders[1].value) do t
    (t1, t)
end

lines!(ax, sols; idxs=1, tspan)
# Drag the slider and errors pop out!

Do you think it's possible to switch back to the more stable APIs for now and switch to the newer experimental Makie APIs once it becomes mature? @asinghvi17

---

Another issue related to this recipe but maybe not to Makie itself: following the way idxs is specified in Plots recipes,

fig, ax, pl = lines(sol, idxs=(1,2,3))

is supposed to plot a 3D figure. However, now in Makie it only shows a 2D figure. If I want a 3D figure, I would need to do

f = Figure()
ax = Axis3(f[1, 1],
   aspect = :data,
)

lines!(ax, sol; idxs=(1,2,3))