Seven Lines of Julia (examples sought)

[JeffreySarnoff]

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.

[JeffreySarnoff]

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.

[rafael.guerra]

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

[mthelm85]

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.

[FedericoStra]

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:

[mforets]

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₂")

@JeffreySarnoff

make a compelling case for Julia to skilled practitioners of whatever is done.

Further details of what is done can be found in ARCH-COMP20 Category Report: Continuous and Hybrid Systems with Nonlinear Dynamics.

[Raf]

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.

                             █ █  ▄▀▀▄             ▄█▄▀█▄▀▀▄      
      ▀▀▀                     ▀    ▀▀             ▀   ▄  ██▀      
    █     █                                         ▄   ██▄       
    ▀     ▀                           ▄▄▄             ▀  ▀     ▄▀▀
      ▀▀▀                                              ▄        ▀▀
                             ██     ██                 █          
                                       ▄           ▄▄▄   ▄▄▄      
                                       █               ▄       ▄▀▄
                                           ██          █       ▀▄▀
   ▄▄                              ██                             
   ▀▀                       ██                      ██   ██       
                                       ▄▄                         
                      ▄▀▄              ▀▀                        ▄
                     ▀▄▄▀       ▄▄▄   ▀▀▄▀ █▄                    ▀
                               █▀▀█ █   ▄▄ ▀                █▄█   
                            ▄   ▀    ▄█▀▄▄▄▄    ██           ▄▄█  
      ▄▄      ▀▀▀▄        ▄█▀▀▄     ▀▄▄▀▄▀  ▄▀█▄▄▄▀         ▄█ ▄ ▀
      ▀▀    █     ▀▄     ▄   ▄▀█    █▄▀ ▀█▄▄    ▀▀           ▀▄▀▄▀
            ▀▄     █   ▄    ▄█▀    ▀ █   ▀▄ █                     
              ▀▄▄▄   ▀▀█   ▀█▄    ▀  ▀▄█   ▀         █▀▀ ▄        
     ▄▀ ▀▄▄              ▄   █▀▄   █ ▀ ▄ ▄█          ▀ ▀▀  ▄▄▄    
  ▄▀▄█▀ ▄█ █              ▄▄  █▀    ██ █▄▀    ▄▀▄        ▄        
 ▀▄ ▄█    ▀                ▀█▀        ▀       ▀▄▀        █        
    ▀▀                            ██                              
       ▄                                                          
      ▀▄█    ▄▄                                                   
             ▀▀                                    ▄              
                      ▀▄                          █ ▀▄            

[Glen]

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

[viralbshah]

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

[juliohm]

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

[rafael.guerra]

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

[lmiq]

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.

[mschauer]

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

[sijo]

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:

[mcreel]

] activate .
instantiate

This is some of my favorite Julia code.

[rafael.guerra]

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

[HenriDeh]

I saw this one on Twitter:

[HenriDeh]

∘ 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.

See Essentials · The Julia Language

[bradeneliason]

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:

[bradeneliason]

I had to submit another one. I’m still leaning about the ApproxFun.jl package, so I’m happy to hear your feedback. This my first real trial of this package other than running a few of the examples.

The following 7 lines (excluding imports) calculates deflection, angle, moment, and shear of a uniformly loaded cantilevered beam and plots the results. This uses ApproxFun.jl find to solution of an ODE. The solution to this ODE is a function which defines the vertical displacement of the beam (under certain assumptions). Figure of the setup given below.

The line beginning B = ... defines the boundary conditions. It’s a fourth order ODE so there are four boundary conditions. The vertical displacement and first derivative (ie, angle) are both zero at the fixed end of the beam (leftendpoint). The moment and shear in the beam (2nd and 3rd) derivatives are zero at the free end of the beam (rightendpoint).

The solution is generated on the next line, v = .... The boundary conditions, B, are all set to 0 using zeros(4).... The differential equation is defined E*I*Derivative()^4 and gets set to a uniform load one(z) (shown as q in the image above).

The rest is just plotting. It’s ugly again to fit into seven lines.

using ApproxFun, Plots
L, E, I = 12.0, 1.0, 1.0
d = 0..L
z = Fun(identity, d)
B = [[Evaluation(d,leftendpoint,k) for k=0:1]... ; [Evaluation(d,rightendpoint,k) for k=2:3]... ;]
v = [B; E*I*Derivative()^4] \ [ zeros(4)..., one(z)]
func_name = zip([v, v', v'', v'''], ["Deflection", "Angle", "Momement", "Shear"])
plot([plot(z, f, title=n, label="") for (f,n) in func_name]..., lw=3)

Result:

It’s easy enough to change the load function:

v = [B; E*I*Derivative()^4] \ [ zeros(4)..., sin(z)]

I haven’t really messed with the units to ensure the solution is scaled correctly, but it does seem to match the analytical solution of the uniformly loaded beam quite well. The comparison:

plot(z, v, lw=3, label="ApproxFun", title="Deflection | Comparison", legend=:left)
plot!(z-> (z^2)*(z^2 - 4*L*z + 6L^2)/(24*E*I), 0, L, lw=3, ls=:dash, label="Analytical")