I am looking for examples: show me some real computational power in seven or fewer lines of Julia that read clearly. It is ok to be using
other packages, though no more than four (and it is ok to use fewer). Each example is intended to make a compelling case for Julia to skilled practitioners of whatever is done.
I was thinking less constructively β what I meant to convey:
There are many expertly crafted packages, a few of which may be essential to demonstrating something important by virtue of their reason for being. To have more than 4 packages used would make it less easy for others to see the underlying powerful simplicity (here of composition as otherβs efforts are working together with oneβs current intent). I appreciate your perspective, and indeed that is something to be shown β¦ perhaps better seen after showing some that require less of the viewer/reader as domain experts are not always equally expert at programming.
Julia code for the Lorenz problem provided in Wikipedia that uses DifferentialEquations.jl.
Lorenzβs equations published in 1963 describe atmospheric convection and are very sensitive to the initial parameters. The work went unnoticed until 1972, when Lorenz gave a conference entitled βPredictability: does the flap of a butterflyβs wings in Brazil set off a tornado in Texas?β, which made famous the βbutterfly effectβ and led to the chaos theory.
The electronic computer used (Royal McBee LGP-30) was according to Lorenz only 1000 times faster than manual computingβ¦ Julia code allows writing this problem compactly with clear syntax, and to collect the solution on todayβs laptops within 1 ms.
using DifferentialEquations, ParameterizedFunctions, Plots
lorenz = @ode_def begin # define the system
dx = Ο*(y - x); dy = x*(Ο - z) - y; dz = x*y - Ξ²*z
end Ο Ο Ξ²
u0=[1.,1.,1.]; tspan=(0.,60.); p=[10.,28.,8/3]; # initial conds., timespan and parms.
prob = ODEProblem(lorenz, u0, tspan, p) # define the problem
sol0 = solve(prob, reltol=1e-9) # solve it
plot(sol0, vars = (1, 2, 3)) # display it
This probably violates the read clearly
part, but:
using CurrentPopulationSurvey
using DataFrames
using DataFramesMeta
oct20 = cpsdata(2020,10) |> DataFrame
unemployment_rate = (@linq oct20 |>
where(map(x -> x in 3:4, :pemlr)) |>
combine(:pwcmpwgt => x -> sum(x) / 10_000))[1,1] /
(@linq oct20 |>
where(map(x -> x in 1:4, :pemlr)) |>
combine(:pwcmpwgt => x -> sum(x) / 10_000))[1,1]
0.06593909683981009
With that, I downloaded the October '20 Current Population Survey public use microdata .zip file from the Census Bureau, unzipped it, parsed it, piped it into a DataFrame, and then computed the non-seasonally-adjusted unemployment rate for the U.S. (go ahead and check it at the Bureau of Labor Statisticsβ data hub ). I think thatβs kind of cool.
The code is a bit golfed and Iβm not exactly sure how to count the lines.
using Plots
function mandelbrot(z) w = z
for n in 1:74
abs2(w) < 4 ? w = w^2 + z : return n
end; 75
end
x, y = range(-0.65, -0.45; length=1600), range(0.51, 0.71; length=1600)
heatmap(x, y, -log.(mandelbrot.(x' .+ y .* im)); aspect_ratio=1)
This is the output:
A visual proof of positive invariance in 7 LOC:
using ReachabilityAnalysis, Plots
@taylorize function vanderpol!(dx, x, p, t)
dx[1] = x[2]; dx[2] = x[2] * (1 - x[1]^2) - x[1]
end
X0 = (1.25 .. 1.55) Γ (2.35 .. 2.45)
prob = @ivp(x' = vanderpol!(x), dim=2, x(0) β X0)
sol = solve(prob, tspan=(0.0, 7.0))
plot(sol, vars=(1, 2), xlabel="xβ", ylabel="xβ")
Further details of what is done can be found in ARCH-COMP20 Category Report: Continuous and Hybrid Systems with Nonlinear Dynamics.
A complete definition of the game of life, and a simulation in the REPL:
using DynamicGrids
const born_survive = (0, 0, 0, 1, 0, 0, 0, 0, 0), (0, 0, 1, 1, 0, 0, 0, 0, 0)
rule = Neighbors(Moore(1)) do neighborhood, state
born_survive[state + 1][sum(neighborhood) + 1]
end
output = REPLOutput(rand(Bool, 100, 100); tspan=1:1000, fps=50, color=:blue)
sim!(output, rule)
(A description of what is going on as requested by @viralbshah)
This cellular automata rule is run over the grid defined by rand(Bool, 100, 100)
, for every cell at every timestep in tstep
. The neighborhood
is an iterable over the surrounding cells. In the most well-known game of life rule, a cell becomes active when the surrounding 8 cells (a βMooreβ neighborhood of radius 1) sum to 2. If the cell is already active it will stay active when the sum is 2 or 3. Swapping the 1s and 0s in born_survive
will give you different GOL rules.
β β ββββ βββββββββ
βββ β ββ β β βββ
β β β βββ
β β βββ β β βββ
βββ β ββ
ββ ββ β
β βββ βββ
β β βββ
ββ β βββ
ββ ββ
ββ ββ ββ ββ
ββ
βββ ββ β
ββββ βββ ββββ ββ β
ββββ β ββ β βββ
β β βββββββ ββ βββ
ββ ββββ βββββ ββββββ βββββββ ββ β β
ββ β ββ β βββ βββ ββββ ββ βββββ
ββ β β βββ β β ββ β
ββββ βββ βββ β βββ β βββ β
ββ βββ β βββ β β β ββ β ββ βββ
βββββ ββ β ββ ββ ββ βββ βββ β
ββ ββ β βββ β βββ β
ββ ββ
β
βββ ββ
ββ β
ββ β ββ
using Measurements, Unitful, Unitful.DefaultSymbols
# Calculate the thickness of tinfoil including uncertainties
Ο = (2.7 Β± 0.2)g/cm^3 # density of aluminum
# For a roll of tinfoil
mass = (250 Β± 10)g
width = (30.5 Β± 0.2)cm
length = (14.24 Β± 0.2)m
# Ο = mass/(width*length*thickness)
@show thickness = ΞΌm(mass/(Ο*width*length))
thickness = ΞΌm(mass / (Ο * width * length)) = 21.3 Β± 1.8 ΞΌm
Iβm loving all the examples being presented here. I would like to encourage folks to submit as many different kinds of examples as possible, and Iβll work on collecting them and getting them onto julialang.org and generally on to social media with something like a #WhyJulia or #WorldOfJulia tag. I think @jiahao used to maintain such a thing a while ago for presentations and such.
It would greatly help to have a 3-5 line description of what it is doing, and why it is interesting, to go along with the code snippet.
-viral
More than 7 lines of code, but still a quite compact pipeline for geospatial machine learning:
using GeoStats, MLJ, CSV, Plots
model = @load DecisionTreeClassifier
df = CSV.read("agriculture.csv") # link to file: https://github.com/JuliaEarth/GeoStats.jl/blob/master/docs/src/data/agriculture.csv
df.crop = categorical(df.crop)
Ξ© = georef(df, (:x,:y))
Ξ©s, Ξ©t = split(Ξ©, 0.2, (1.0,-1.0))
π― = ClassificationTask([:band1,:band2,:band3,:band4], :crop)
π« = LearningProblem(Ξ©s, Ξ©t, π―)
π = PointwiseLearn(model)
Ξ©Μt = solve(π«, π)
plot(Ξ©Μt, ms=0.2, c=:cividis, title="crop (prediction)")
Commented example is here: https://juliaearth.github.io/GeoStats.jl/stable/workflow/#Defining-problems
Like with Julia, experts say that Shakespearean language can be tricky but that acting it out makes it fun and understandable. On the same vein, the following MIT list of Shakespearean insults might be a good way to learn both languages:
using CSV, DataFrames, HTTP
f = CSV.File(HTTP.get("http://web.mit.edu/dryfoo/Funny-pages/shakespeare-insult-kit.html").body);
sk = DataFrame(f[17:66]) # TODO: findfirst "artless" word to set range
r = size(sk,1)
sk = DataFrame([split(sk[i,1]) for i in 1:r]) # output = 3 rows x 50 columns
Shakespeare() = println("Thou "*sk[1,rand(1:r)]*" "*sk[2,rand(1:r)]*" "*sk[3,rand(1:r)])
Then call Shakespeare:
julia> Shakespeare()
Thou artless beef-witted joithead
julia> Shakespeare()
Thou villainous motley-minded giglet
With 15 lines I could write a code to perform a particle simulation with periodic boundary conditions, a langevin thermostat, and a quadratic potential between the particles, and produce an animation:
using Plots ; ENV["GKSwstype"]="nul"
const N, Ο, Ξt, Ξ», T, k = 100, 1000, 0.01, 1e-3, 0.26, 1e-6
const x, v, f = -0.5 .+ rand(3,N), -0.01 .+ 0.02*randn(3,N), zeros(3,N)
wrap(x,y) = (x-y) > 0.5 ? (x-y)-1 : ( (x-y) < -0.5 ? (x-y)+1 : (x-y) )
anim = @animate for t in 1:Ο
f .= 0
for i in 1:N-1, j in i+1:N
f[:,i] .+= wrap.(x[:,i],x[:,j]) .- Ξ» .* v[:,i]
f[:,j] .+= wrap.(x[:,j],x[:,i]) .- Ξ» .* v[:,j]
end
x .= wrap.(x .+ v*Ξt .+ (f/2)*Ξt^2,zeros(3))
v .= v .+ f*Ξt .+ sqrt.(2*Ξ»*k*T*Ξt)*randn()
scatter(x[1,:],x[2,:],x[3,:],label="",lims=(-0.5,0.5),aspect_ratio=1,framestyle=:box)
end
gif(anim,"anim.gif",fps=10)
As Tamas mentioned, if we go to a few more lines, this could include an actual Lennard-Jones potential and be parallelized. And that without using any package besides Plots.
Infinizoom! (Background: Simulates a Brownian motion and uses the Brownian self-similarity Wiener process - Wikipedia combined with simple data-augmentation/Brownian bridge sampling to zoom in infinitely.)
using Makie
dt = 0.0005; Ο΅ = 0.1; t = collect(-1.0:dt:1.0); n = length(t)
x = sqrt(dt)*cumsum(randn(n)); x .-= x[endΓ·2]
T, X = Node(t), Node(x)
p = lines(T, X); display(p)
xlims!(-1.0, 1.0); ylims!(-1.5,1.5)
while true
t .*= (1 + Ο΅)
x .*= sqrt(1 + Ο΅)
if t[end] > 2.0
t .= collect(-1.0:dt:1.0)
x[1:2:end] = x[nΓ·4+1:3nΓ·4+1]
x[2:2:end] = (x[1:2:end-2]+x[3:2:end])/2 + sqrt(dt)/2*randn(nΓ·2)
x[endΓ·2+1] = 0
end
T[] = t; X[] = x; sleep(0.02)
end
Nice! I was thinking of doing something similar, but I couldnβt find the time, nor fit it in less than 7 lines (you cheated! ).
I think that, to make it better, the number of iterations in mandelbrot
has to increase the more one zooms in, otherwise the boundary of the fractal becomes too imprecise. Here is a complete version that does so:
using Plots
exprange(start, stop, len) = exp.(range(log(start), log(stop), length=len))
function mandelbrot(z; lim=75) w = z
for n = 1:lim; abs2(w) < 4 ? w = w^2 + z : return n end
lim + 1
end
xβ, yβ = -0.5626805, 0.6422555
anim = @animate for (r, l) in zip(exprange(2, 1.35e-6, 120), exprange(100, 2500, 120))
x = range(xβ-r, xβ+r; length=600); y = range(yβ-r, yβ+r; length=600);
heatmap(x, y, -log.(log.(mandelbrot.(x' .+ y .* im; lim=round(l))));
legend=:none, border=:none, ticks=:none, size=(600,600), ratio=1)
end
g = gif(anim; fps=12)
Apparently I cannot attach the animated GIF file (which you can find here), so I include only the last frame:
Calculating a cross-correlation in any number of dimensions:
for I in CartesianIndices(img)
for J in CartesianIndices(kernel) # kernel with centered indices
if I+J in CartesianIndices(img)
filtered[I] += img[I+J] * kernel[J]
end
end
end
(Inspired by Multidimensional algorithms and iteration and Knowing where you are: custom array indices in Julia)
The code in action:
using OffsetArrays
using Makie
# A 3D "image" of random values (either 0 or 1)
img = rand([0.0, 1.0], 20, 60, 40)
# Plot
s1 = volume(img, algorithm=:iso, resolution=(800, 500))
# Function to generate an N-dimensional kernel of uniform values that sum to 1
# The center of the hypercube will have indices (0, 0, ..., 0)
uniform_kernel(n, dim) = OffsetArray(fill(1/n^dim, fill(n, dim)...), fill(-nΓ·2-1, dim)...)
# 5x5x5 array of uniform values that sum to 1
kernel = uniform_kernel(5, 3) # 3 dimensions
filtered = zero(img)
for I in CartesianIndices(img)
for J in CartesianIndices(kernel)
if I+J in CartesianIndices(img)
filtered[I] += img[I+J] * kernel[J]
end
end
end
s2 = volume(filtered, algorithm=:iso, resolution=(800, 500))
Original image:
Filtered image:
] activate .
instantiate
This is some of my favorite Julia code.
MΓΆbius discovered the one-sided non-orientable strip in 1858 while serving as the chair of astronomy at the University of Leipzig. The concept inspired mathematicians in the development of the field of topology as well as artists like M.C. Escher.
Fifteen lines of Julia were required to compute and plot the Mobius strip and a moving normal vector, which changes color after each 2Ο
revolution. This exceeds the limit but someone might be able to compress it further.
(PS: edited code for better readability. Without sleep
, the last for
loop could be written as a comprehension too)
using Makie, ForwardDiff, LinearAlgebra
const FDd = ForwardDiff.derivative
M(u,v) = ((r + v/2*cos(u/2))*cos(u), (r + v/2*cos(u/2))*sin(u), v/2*sin(u/2))
βMβu(u,v) = [FDd(u->M(u,v)[1],u), FDd(u->M(u,v)[2],u), FDd(u->M(u,v)[3],u)]
βMβv(u,v) = [FDd(v->M(u,v)[1],v), FDd(v->M(u,v)[2],v), FDd(v->M(u,v)[3],v)]
r = 1.0; w = 0.5; m = 90; clr = (:cyan,:red,:cyan,:red,:cyan)
u = range(0, 8Ο, length = 4m); v = range(-w, w, length = 3)
N = [ -cross(βMβu(ui, 0.0), βMβv(ui, 0.0))[j] for ui in u, j in 1:3 ]
N = N ./ sqrt.(sum(abs2,N,dims=2)); N = vec.(Point3f0.(N[:,1],N[:,2],N[:,3]))
xs, ys, zs = [[p[i] for p in M.(u, v')] for i in 1:3]
P0 = vec.(Point3f0.(xs[:,2], ys[:,2], zs[:,2]))
surface(xs[1:m+1,:], ys[1:m+1,:], zs[1:m+1,:], backgroundcolor=:black)
for i in 1:4m
sleep(0.01); k = 1 + iΓ·(m+0.5);
arrows!([P0[i]],[N[i]],lengthscale=0.3,arrowsize=0.05,arrowcolor=clr[k],linecolor=clr[k])
end
β is the function composition function. β(f,g,h)(x)
turns into f(g(h(x)))
.
fill(F,n)
populates an array with n times the function F. They are composed by β to form the function
x -> F(F(F(F...F(x)))...)
.
Thatβs is, you just call it with any argument you like.
Hereβs my seven lines (excluding imports). It calculates and plots pursuit curves using DifferentialEquations.jl. Iβve renamed the pursuer and pursued to fox
and rabbit
for clarity. Itβs kind of fun to play around with different paths for the rabbit to take. I have a parameter, k
, which alters the speed of the fox.
The fourth line beginning Rx, Ry,...
is a little goofy to fit within the space constraint, but it works and makes what Iβm plotting clear as day.
using Plots, DifferentialEquations, LinearAlgebra
rabbit(t) = [cos(t), sin(t)] # Rabbit Running in a Circle
fox(u, k, t) = k * (rabbit(t) - u) / norm(rabbit(t)-u) # Fox chases rabbit at speed k
prob = ODEProblem(fox, [2.0, 0.0], (0.0, 10), 0.8)
sol = solve(prob, saveat=0.1);
Rx, Ry, Fx, Fy = [getindex.(vcat.(rabbit.(sol.t), sol.u), i) for i=1:4]
plot(Rx, Ry, label="Rabbit", lw=sol.t)
plot!(Fx, Fy, label="Fox", lw=sol.t, aspect_ratio=:equal)
Hereβs a bonus animation: