Making animated and interactive scientific visualizations in Makie.jl

makie_tutorial.jl

# # Making animated and interactive scientific 
# # visualizations in Makie.jl

# Alrighty, in this short tutorial I'll teach you how to make
# animated scientific such as this video, or this application.

# Prerequisitives:
# * Basic familiarity with Julia (obviously)
# * Basic familiarity with Makie.jl

# # The process
# 0. Learn about `Observable`s
# 1. Initialize simulation in a stepping manner
# 2. Initialize the `Observable`s of the animation 
# 3. Plot the `Observable`s and any other static elements
# 4. Create the "animation stepping function"
# 5. Test it
# 6. Save animations to videos
# 7. Interactive application

# Let's start!
# Load packages that will be used
using DynamicalSystems
using OrdinaryDiffEq
using GLMakie
using DataStructures: CircularBuffer

# %% 0. Learn about `Observable`s
# An `Observable` is a mutable container of an object of type `T`.
# `T` can be any type. The value of an `Observable` can then be 
# changed on the spot, just like updating any mutable container.
# (This is similar to the base structure `Ref`, if you're familiar)

# The important part here is that `Observable`s can be "listened" to.
# What does this mean...?

o = Observable(1) # Observable with values of type `Int`

l1 = on(o) do val # Create a listener `l1` of observable.
    println("Observable now has value $val")
end

# `l1` is triggered each time the value of `o` is updated.
# (demo in REPL, set `o[] = 2`.)


# We care about `Observable`s because Makie.jl is hooked up
# to this "listener" system. If any plotted element is
# initialized as an observable, then Makie.jl understands this.
# Updating the observable would update the plot values.

# For example:
ox = 1:4
oy = Observable(rand(4))
lw = Observable(2)

fig, ax = lines(ox, oy; linewidth = lw)
ylims!(ax, 0, 1)

lw[] = 50
oy[] = rand(4)

# This simple process is the basis of creating 
# both animations, as well as interactive apps with Makie.jl.
# So in the following let's begin making a simple visualization
# and interactive application of the double pendulum


# %% 1. Initialize simulation in a stepping manner
# (this can also be done with a pre-run simulation)
# the goal is to create a "step" function which
# once called it progresses the data for one animation frame
const L1 = 1.0
const L2 = 0.9
M = 2
u0 = [π/3, 0, 3π/4, -2]
dp = Systems.double_pendulum(u0; L1, L2)

# Solve diffeq with constant step for smoother curves
diffeq = (alg = Tsit5(), adaptive = false, dt = 0.005)

integ = dp.integ

function xycoords(state)
    θ1 = state[1]
    θ2 = state[3]
    x1 = L1 * sin(θ1)
    y1 = -L1 * cos(θ1)
    x2 = x1 + L2 * sin(θ2)
    y2 = y1 - L2 * cos(θ2)
    return x1,x2,y1,y2
end

# `integ` is an integrator. `step!(integ)` progresses the integrator
# for one step. `integ.u` is the system state at current step. 
# Then `xycoords` converts the states of the integrator
# to their plottable format. So we can imagine something like 
function progress_for_one_step!(integ)
    step!(integ)
    return xycoords(integ)
end
# to be our stepping function that returns the new data.

# If we had finite data instead of a forever-running animation, 
# then the "stepping function" would simply be to progress the index `i`
# of the existing data one step forwards...


# %% 2. Initialize the `Observable`s of the animation 
# You need to think of this in advance: what things will to be 
# animated, and what other plotting elements will be static? 
# Animated elements will need to become `Observable`s.

# Here the animated elements will be: balls and rods making the
# double pendulum, and the tail (trajectory) of the pendulum.
x1,x2,y1,y2 = xycoords(u0)
rod   = Observable([Point2f(0, 0), Point2f(x1, y1), Point2f(x2, y2)])
balls = Observable([Point2f(x1, y1), Point2f(x2, y2)])
# (Remember: the most optimal way to plot 2D things in Makie.jl is to
# give it a vector of `Point2f`, the coordinates for the plot)

# Here we have initialized two _different_ observables, because
# rods and balls will be plotted in a different manner (lines/scatter)

# Next is the observable for the tail
tail = 300 # length of plotted trajectory, in units of `dt`
# The circular buffer datastructure makes making stepping-based
# animations very intuitive
traj = CircularBuffer{Point2f}(tail)
fill!(traj, Point2f(x2, y2)) # add correct values to the circular buffer
traj = Observable(traj) # make it an observable


# %% 3. Plot the `Observable`s and any other static elements
# Before plotting we need to initialie a figure
fig = Figure(); display(fig)
# in my experience it leads to cleaner code if we first initialize 
# an axis and populate it accordingly.
ax = Axis(fig[1,1])

# Now we plot the observables _directly_! First the pendulum
lines!(ax, rod; linewidth = 4, color = :purple)
scatter!(ax, balls; marker = :circle, strokewidth = 2, 
    strokecolor = :purple,
    color = :black, markersize = [8, 12]
)

# then its trajectory, with a nice fadeout color
c = to_color(:purple)
tailcol = [RGBAf(c.r, c.g, c.b, (i/tail)^2) for i in 1:tail]
lines!(ax, traj; linewidth = 3, color = tailcol)

# We can also plot now any other static elements
ax.title = "double pendulum"
ax.aspect = DataAspect()
l = 1.05(L1+L2)
xlims!(ax, -l, l)
ylims!(ax, -l, 0.5l)

# %% 4. Create the "animation stepping function"
# Using the functions of step 1, we now define a function
# that updates the observables. Makie.jl understands observable
# updates and directly reflects this on the plotted elements.
function animstep!(integ, rod, balls, traj)
    x1,x2,y1,y2 = progress_for_one_step!(integ)
    rod[] = [Point2f(0, 0), Point2f(x1, y1), Point2f(x2, y2)]
    balls[] = [Point2f(x1, y1), Point2f(x2, y2)]
    push!(traj[], Point2f(x2, y2))
    traj[] = traj[] # <- important! Updating in-place the value of an
                    # `Observable` does not trigger an update!
end

# %% 5. Test it
for i in 1:1000
    animstep!(integ, rod, balls, traj)
    sleep(0.001)
end

# cool it works. Let's wrap up the creation of the observables
# and plots in a function (just to re-initialie everything)
function makefig(u0)
    dp = Systems.double_pendulum(u0; L1, L2)
    integ = dp.integ
    x1,x2,y1,y2 = xycoords(u0)
    rod   = Observable([Point2f(0, 0), Point2f(x1, y1), Point2f(x2, y2)])
    balls = Observable([Point2f(x1, y1), Point2f(x2, y2)])
    traj = CircularBuffer{Point2f}(tail)
    fill!(traj, Point2f(x2, y2)) # add correct values to the circular buffer
    traj = Observable(traj) # make it an observable
    fig = Figure(); display(fig)
    ax = Axis(fig[1,1])
    lines!(ax, rod; linewidth = 4, color = :purple)
    scatter!(ax, balls; marker = :circle, strokewidth = 2, 
        strokecolor = :purple,
        color = :black, markersize = [8, 12]
    )
    lines!(ax, traj; linewidth = 3, color = tailcol)
    ax.title = "double pendulum"
    ax.aspect = DataAspect()
    l = 1.05(L1+L2)
    xlims!(ax, -l, l)
    ylims!(ax, -l, 0.5l)
    # also return the figure object, we'll ues it!
    return fig, integ, rod, balls, traj
end
    
# %% 6. Save animations to videos
fig, integ, rod, balls, traj = makefig(u0)
frames = 1:200
record(fig, "video.mp4", frames; framerate = 60) do i # i = frame number
    for j in 1:5 # step 5 times per frame
        animstep!(integ, rod, balls, traj)
    end
    # any other manipulation of the figure here...
end # for each step of this loop, a frame is recorded

# %% 7. Interactive application
# Makie.jl has tremendously strong capabilities for real-time
# interactivity. To learn all of this takes time of course,
# and you'll need to consult the online documentation.
# Here we will do two interactions: 1) a play/stop button
# 2) clicking on the screen and getting a new initial condition!

fig, integ, rod, balls, traj = makefig(u0)
# The run button is actually pretty simple, we'll add it below the plot
run = Button(fig[2,1]; label = "run", tellwidth = false)
# This button will start/stop an animation. It's actually surprisingly
# simple to do this. The magic code is:
isrunning = Observable(false)
on(run.clicks) do clicks; isrunning[] = !isrunning[]; end
on(run.clicks) do clicks
    @async while isrunning[]
        isopen(fig.scene) || break # ensures computations stop if closed window
        animstep!(integ, rod, balls, traj)
        sleep(0.001) # or `yield()` instead
    end
end

# `on` an important Observables function when one starts
# doing advanced stuff. It triggers a piece of code once an observable
# triggers its update.

# We'll add one more interactive feature which will trigger once
# we click on the axis. Notice that by default makie performs a zoom
# once one clicks on the axis, so we'll disable this
ax = content(fig[1,1])
Makie.deactivate_interaction!(ax, :rectanglezoom)
# and we'll add a new trigger using the `select_point` function:
spoint = select_point(ax.scene, marker = :circle)

# Now we see that we can click on the screen and the `spoint` updates!

# Okay, let's do something useful when it triggers
function θωcoords(x, y)
    θ = atan(y,x) + π/2
    return SVector(θ,0,0,0)
end

on(spoint) do z
    x, y = z
    u = θωcoords(x, y)
    reinit!(integ, u)
    # Reset tail and balls to new coordinates
    x1,x2,y1,y2 = xycoords(u)
    traj[] .= fill(Point2f(x2, y2), length(traj[]))
    traj[] = traj[]
    rod[] = [Point2f(0, 0), Point2f(x1, y1), Point2f(x2, y2)]
    balls[] = [Point2f(x1, y1), Point2f(x2, y2)]
end

Nice demo! 😄

@Datseris Point2f0 and RGBAf0 have been renamed to Point2f and RGBAf, in case you want to update the script.

thanks, just updated!

I don't know how to make a pull request of a gist, but here's an updated code that works:

# # Making animated and interactive scientific 
# # visualizations in Makie.jl

# Alrighty, in this short tutorial I'll teach you how to make
# animated scientific such as this video, or this application.

# Prerequisitives:
# * Basic familiarity with Julia (obviously)
# * Basic familiarity with Makie.jl

# # The process
# 0. Learn about `Observable`s
# 1. Initialize simulation in a stepping manner
# 2. Initialize the `Observable`s of the animation 
# 3. Plot the `Observable`s and any other static elements
# 4. Create the "animation stepping function"
# 5. Test it
# 6. Save animations to videos
# 7. Interactive application

# Let's start!
# Load packages that will be used
using DynamicalSystems
using OrdinaryDiffEq
using GLMakie
using DataStructures: CircularBuffer

# %% 0. Learn about `Observable`s
# An `Observable` is a mutable container of an object of type `T`.
# `T` can be any type. The value of an `Observable` can then be 
# changed on the spot, just like updating any mutable container.
# (This is similar to the base structure `Ref`, if you're familiar)

# The important part here is that `Observable`s can be "listened" to.
# What does this mean...?

o = Observable(1) # Observable with values of type `Int`

l1 = on(o) do val # Create a listener `l1` of observable.
    println("Observable now has value $val")
end

# `l1` is triggered each time the value of `o` is updated.
# (demo in REPL, set `o[] = 2`.)


# We care about `Observable`s because Makie.jl is hooked up
# to this "listener" system. If any plotted element is
# initialized as an observable, then Makie.jl understands this.
# Updating the observable would update the plot values.

# For example:
ox = 1:4
oy = Observable(rand(4))
lw = Observable(2)

fig, ax = lines(ox, oy; linewidth = lw)
ylims!(ax, 0, 1)

lw[] = 50
oy[] = rand(4)

# This simple process is the basis of creating 
# both animations, as well as interactive apps with Makie.jl.
# So in the following let's begin making a simple visualization
# and interactive application of the double pendulum


# %% 1. Initialize simulation in a stepping manner
# (this can also be done with a pre-run simulation)
# the goal is to create a "step" function which
# once called it progresses the data for one animation frame
const L1 = 1.0
const L2 = 0.9
M = 2
u0 = [π/3, 0, 3π/4, -2]
dp = Systems.double_pendulum(u0; L1, L2)

# Solve diffeq with constant step for smoother curves
diffeq = (alg = Tsit5(), adaptive = false, dt = 0.005)

integ = dp.integ

function xycoords(state)
    θ1 = state[1]
    θ2 = state[3]
    x1 = L1 * sin(θ1)
    y1 = -L1 * cos(θ1)
    x2 = x1 + L2 * sin(θ2)
    y2 = y1 - L2 * cos(θ2)
    return x1,x2,y1,y2
end

# `integ` is an integrator. `step!(integ)` progresses the integrator
# for one step. `integ.u` is the system state at current step. 
# Then `xycoords` converts the states of the integrator
# to their plottable format. So we can imagine something like 
function progress_for_one_step!(integ)
    step!(integ)
    return xycoords(integ)
end
# to be our stepping function that returns the new data.

# If we had finite data instead of a forever-running animation, 
# then the "stepping function" would simply be to progress the index `i`
# of the existing data one step forwards...


# %% 2. Initialize the `Observable`s of the animation 
# You need to think of this in advance: what things will to be 
# animated, and what other plotting elements will be static? 
# Animated elements will need to become `Observable`s.

# Here the animated elements will be: balls and rods making the
# double pendulum, and the tail (trajectory) of the pendulum.
x1,x2,y1,y2 = xycoords(u0)
rod   = Observable([Point2f(0, 0), Point2f(x1, y1), Point2f(x2, y2)])
balls = Observable([Point2f(x1, y1), Point2f(x2, y2)])
# (Remember: the most optimal way to plot 2D things in Makie.jl is to
# give it a vector of `Point2f`, the coordinates for the plot)

# Here we have initialized two _different_ observables, because
# rods and balls will be plotted in a different manner (lines/scatter)

# Next is the observable for the tail
tail = 300 # length of plotted trajectory, in units of `dt`
# The circular buffer datastructure makes making stepping-based
# animations very intuitive
traj = CircularBuffer{Point2f}(tail)
fill!(traj, Point2f(x2, y2)) # add correct values to the circular buffer
traj = Observable(traj) # make it an observable


# %% 3. Plot the `Observable`s and any other static elements
# Before plotting we need to initialie a figure
fig = Figure(); display(fig)
# in my experience it leads to cleaner code if we first initialize 
# an axis and populate it accordingly.
ax = Axis(fig[1,1])

# Now we plot the observables _directly_! First the pendulum
lines!(ax, rod; linewidth = 4, color = :purple)
scatter!(ax, balls; marker = :circle, strokewidth = 2, 
    strokecolor = :purple,
    color = :black, markersize = [8, 12]
)

# then its trajectory, with a nice fadeout color
c = to_color(:purple)
tailcol = [RGBAf(c.r, c.g, c.b, (i/tail)^2) for i in 1:tail]
lines!(ax, traj; linewidth = 3, color = tailcol)

# We can also plot now any other static elements
ax.title = "double pendulum"
ax.aspect = DataAspect()
l = 1.05(L1+L2)
xlims!(ax, -l, l)
ylims!(ax, -l, 0.5l)

# %% 4. Create the "animation stepping function"
# Using the functions of step 1, we now define a function
# that updates the observables. Makie.jl understands observable
# updates and directly reflects this on the plotted elements.
function animstep!(integ, rod, balls, traj)
    x1,x2,y1,y2 = progress_for_one_step!(integ)
    rod[] = [Point2f(0, 0), Point2f(x1, y1), Point2f(x2, y2)]
    balls[] = [Point2f(x1, y1), Point2f(x2, y2)]
    push!(traj[], Point2f(x2, y2))
    traj[] = traj[] # <- important! Updating in-place the value of an
                    # `Observable` does not trigger an update!
end

# %% 5. Test it
for i in 1:1000
    animstep!(integ, rod, balls, traj)
    sleep(0.001)
end

# cool it works. Let's wrap up the creation of the observables
# and plots in a function (just to re-initialie everything)
function makefig(u0)
    dp = Systems.double_pendulum(u0; L1, L2)
    integ = dp.integ
    x1,x2,y1,y2 = xycoords(u0)
    rod   = Observable([Point2f(0, 0), Point2f(x1, y1), Point2f(x2, y2)])
    balls = Observable([Point2f(x1, y1), Point2f(x2, y2)])
    traj = CircularBuffer{Point2f}(tail)
    fill!(traj, Point2f(x2, y2)) # add correct values to the circular buffer
    traj = Observable(traj) # make it an observable
    fig = Figure(); display(fig)
    ax = Axis(fig[1,1])
    lines!(ax, rod; linewidth = 4, color = :purple)
    scatter!(ax, balls; marker = :circle, strokewidth = 2, 
        strokecolor = :purple,
        color = :black, markersize = [8, 12]
    )
    lines!(ax, traj; linewidth = 3, color = tailcol)
    ax.title = "double pendulum"
    ax.aspect = DataAspect()
    l = 1.05(L1+L2)
    xlims!(ax, -l, l)
    ylims!(ax, -l, 0.5l)
    # also return the figure object, we'll ues it!
    return fig, integ, rod, balls, traj
end
    
# %% 6. Save animations to videos
fig, integ, rod, balls, traj = makefig(u0)
frames = 1:200
record(fig, "video.mp4", frames; framerate = 60) do i # i = frame number
    for j in 1:5 # step 5 times per frame
        animstep!(integ, rod, balls, traj)
    end
    # any other manipulation of the figure here...
end # for each step of this loop, a frame is recorded

# %% 7. Interactive application
# Makie.jl has tremendously strong capabilities for real-time
# interactivity. To learn all of this takes time of course,
# and you'll need to consult the online documentation.
# Here we will do two interactions: 1) a play/stop button
# 2) clicking on the screen and getting a new initial condition!

fig, integ, rod, balls, traj = makefig(u0)
# The run button is actually pretty simple, we'll add it below the plot
run = Button(fig[2,1]; label = "run", tellwidth = false)
# This button will start/stop an animation. It's actually surprisingly
# simple to do this. The magic code is:
isrunning = Observable(false)
on(run.clicks) do clicks; isrunning[] = !isrunning[]; end
on(run.clicks) do clicks
    @async while isrunning[]
        isopen(fig.scene) || break # ensures computations stop if closed window
        animstep!(integ, rod, balls, traj)
        sleep(0.001) # or `yield()` instead
    end
end

# `on` an important Observables function when one starts
# doing advanced stuff. It triggers a piece of code once an observable
# triggers its update.

# We'll add one more interactive feature which will trigger once
# we click on the axis. Notice that by default makie performs a zoom
# once one clicks on the axis, so we'll disable this
ax = content(fig[1,1])
Makie.deactivate_interaction!(ax, :rectanglezoom)
# and we'll add a new trigger using the `select_point` function:
spoint = select_point(ax.scene, marker = :circle)

# Now we see that we can click on the screen and the `spoint` updates!

# Okay, let's do something useful when it triggers
function θωcoords(x, y)
    θ = atan(y,x) + π/2
    return SVector(θ,0,0,0)
end

on(spoint) do z
    x, y = z
    u = θωcoords(x, y)
    reinit!(integ, u)
    # Reset tail and balls to new coordinates
    x1,x2,y1,y2 = xycoords(u)
    traj[] .= fill(Point2f(x2, y2), length(traj[]))
    traj[] = traj[]
    rod[] = [Point2f(0, 0), Point2f(x1, y1), Point2f(x2, y2)]
    balls[] = [Point2f(x1, y1), Point2f(x2, y2)]
end```

Also, I built a little DS explorer tool with a bifurcation diagram that you can click and return maps, but it crashes all the time, intermittently enough that I can't make a MWE. Do you have any experience with this? I have no stacktrace or logs to troubleshoot from and I'm at a loss.

thanks updated!

Open an issue at InteractiveDynamics.jl and attach your working code and I'll see if I can have a look.

Regarding the @async loop above, I wonder if there is the possibility that two button clicks in rapid succession could spawn a second while loop without terminating the first one.

As an alternative, how about the code pattern below? It maintains a persistent async while loop. If the simulation is stopped, the loop goes to sleep via wait, and can be reawakened through notify on a Condition.

using GLMakie, Random

spins = Observable(zeros(10, 10))

fig = Figure()
ax = Axis(fig[1,1])

heatmap!(ax, spins)


# Async loop to perform simulations while `isrunning` is true. Need a Condition
# to wake up a thread on GUI events.
isrunning = Observable(false)
isrunning_notifier = Condition()
on(cond -> cond && notify(isrunning_notifier), isrunning)
errormonitor(@async while true
    if isrunning[]
        randn!(spins[])
        notify(spins)
        yield()
    else
        wait(isrunning_notifier)
    end
end)

# Button to start/stop animation
label = map(cond -> cond ? "Stop" : "Run", isrunning)
runbutton = Button(fig[2,1]; label, tellwidth = false)
on(runbutton.clicks) do _
    isrunning[] = !isrunning[]
end

# Stop animation if scene is not visible
on(fig.scene.visible) do visible
    isrunning[] = visible
end

display(fig)

thanks for the suggetsion!

I realized today that it's good to wrap the @async block inside an errormonitor. This way if the simulation crashes, the error message won't be hidden. Example above updated.

https://docs.julialang.org/en/v1/base/parallel/#Base.errormonitor

thanks a lot for these recommendations. When I make an updated tutorial for animations (once finally Makie 1.0 is released) I will include all your suggestions!