[show and tell] Blog post on interval arithmetic in Pluto.jl

I just wrote a blog post entirely in Julia and Pluto.jl! It’s the story of how I flew too close to the sun implementing my own interval arithmetic library, and lessons learned along the way. It’s got some math, some twists, love interests, action. A blog post for the whole family.

This is the first blog post I’ve ever published since my prehistoric wordpress days, and I have to say Pluto.jl was a dream to work with! It’s so good that it’s making me start to consider Julia for projects that I never would’ve otherwise done in Julia or in notebooks for that matter. It doesn’t have any of the same is-this-cell-stale anxiety that I had with Jupyter. Shout out to @fonsp and team!

27 Likes

Nice. Very clean. How did you host it? Is this more than just slapping the generated HTML onto a static site?

Thanks! I just exported to HTML and threw it onto a static site! I don’t know if this is the best solution, but it works for now…

5 Likes

It looks great. I guess I’m just amazed that it’s that easy to make that good looking of a site with Pluto. It makes it look like you did a ton of work, putting a working Binder integration button on the top and everything :sweat_smile:

3 Likes

Nice post!

Let me share another example to further illustrate your statement that:

Practically speaking it means that although interval arithmetic generally works very well without supervision, you may want to write hand-crafted interval arithmetic functions for certain operations or consider different possible rewrites of your computation

This small example mixes intervals, the sine function… and differential equations, for the whole Julia :juliabouncer: familly. (Credits to @AnderGray for this nice example).

The function x(t) = x_0 \exp\left(\cos(t) - 1\right), x_0 \in \mathbb{R}, t \geq 0 solves the differential equation x'(t) = -x(t) ~ \sin(t).

Standard integration schemes fail to produce helpful solutions if the initial state is an interval, for reasons explained in your post, eg. dependency problems. This is easy to check:

using DifferentialEquations, IntervalArithmetic

# initial condition
x₀ = [-1 .. 1]

# define the problem
function f(dx, x, p, t)
    dx[1] = -x[1] * sin(t)
end

# pass to solvers
prob = ODEProblem(f, x₀, (0.0, 2.0))
sol = solve(prob, Tsit5(), adaptive=false, dt=0.05, reltol=1e-6);

# there is no plot recipe readily available so we create a plot with the help of LazySets.jl
using LazySets, Plots
using LazySets: Interval

out = [Interval(sol.t[i]) × Interval(sol.u[i][1]) for i in 1:20]

plot(out, vars=(0, 1), xlab="t", ylab="x(t)", lw=3.0, alpha=1., c=:black, marker=:none, title="Standard integrator with interval initial condition")

See that the solution explodes already for 1 second. On the other hand, specialized algorithms can handle this case without noticeable wrapping effect, producing a sequence of sets whose union covers the true solution for all initial points.

using ReachabilityAnalysis, Plots

# define the model
@taylorize function f(dx, x, p, t)
    dx[1] = -x[1] * sin(t)
end

# define the set of initial states; in this case, an interval
X0 = -1 .. 1

# define the initial-value problem
prob = @ivp(x' = f(x), x(0) ∈ X0, dim=1);

# solve it
sol = solve(prob, alg=TMJets21a(abstol=1e-10), T=15);

# it is illustrative to plot the computed flowpipe and the known analytic solution
# for a range of initial conditions that cover the range of `U0`.
analytic_sol(x0) = t -> x0 * exp(cos(t) - 1.0)

dt = range(0, 15, length=100)
x0vals = range(-1, 1, length=25)

fig = plot()
plot!(fig, sol, vars=(0, 1), lw=0.0, title="Specialized (Taylor-model based) integrator")
[plot!(fig, dt, analytic_sol(x0).(dt), lab="", c=:magenta, lw=2.0, xlab="t", ylab="x(t)") for x0 in x0vals]; #!jl
fig

We will discuss this and other examples at our JuliaCon workshop “It’s all set: a hands-on introduction to JuliaReach” next monday :julia:

7 Likes

Nice! What’s the trick here? What does @taylorize do?

Do you mean the math? You could check works by Berz and Makino, M. Joldes, X. Chen, F. Bünger to name just a few.

About the implementation, you could check a JuliaCon 2018 presentation by @lbenet and @dpsanders about Taylor models. The methods used above to propagate solutions of nonlinear ODEs is more recent though.

It is an optimization strategy used in TaylorIntegration.jl; I checked removing it in this example and there is no difference with and without it.

2 Likes

Great post!

For interactivity, you could use GitHub - JuliaPluto/PlutoSliderServer.jl: Web server to run just the @bind parts of a Pluto.jl notebook in addition. With it, users can manipulate PlutoUI elements on the web page itself (without running Pluto locally) and the shown notebook content updates in real-time.
The disadvantage is that it needs som (virtual) machine to run the server.

PlutoSliderServer is used e.g. here:

https://computationalthinking.mit.edu/Spring21/images/

Thanks! Yeah, I would’ve liked to add some interactivity but unfortunately the site is fully static, so I don’t have a VM to spare. It would really be nice if Julia could run in the browser via WebAssembly or something instead.

Wasm support for Julia would be great!
There are some efforts in this direction: GitHub - Keno/julia-wasm: Running julia on wasm, Web Platform Projects – Summer of Code

2 Likes