Seven Lines of Julia (examples sought)

Very nice :smile: Back in the day when I needed continued fractions, I also added \structs to the TeX code, so that the numbers didn’t keep getting smaller.

Edit: It was very politely pointed out to me that it’s de rigeur in this thread to provide code. Although it’s not only Julia that lets you do this, I really like the swap semantics in Julia. The following is code that computes a random permutation of an input vector in six lines:

function random_permutation(a)
    for m = 1:length(a)-1
        l = m+rand(0:length(a)-m)
        a[l],a[m]=a[m],a[l]
    end
end

As part of a project where I’m porting some FORTRAN (FORTRAN IV!) code to Julia, I really like that Julia is not only much more succinct but also much more general.

1 Like

Dilbert comic strips do not need introduction:

using Gumbo, HTTP, FileIO, Dates, Plots; gr()

function dilbert(t::Date = today() - Day(1))
    r = HTTP.get("https://dilbert.com/strip/"*string(t))
    r_parsed = parsehtml(String(r.body))
    webaddress = r_parsed.root[1][26].attributes["content"]
    download(webaddress) |> load |> x -> display(plot(x,dpi=300, ticks=nothing))
end
dilbert(s::String) = let t = rand(Date(1990,1,1):Day(1):today() - Day(1)); dilbert(t); println(t) end

Then call: dilbert()

Or: dilbert(Date(2021,06,26))

Or better, do repeat calls to: dilbert("rand")

NB: edited random function, thanks to @Syx_Pek

16 Likes

I recommend the following improvement to the last line.

dilbert(s::String) = let t = rand(Date(1990,1,1):Day(1):today() - Day(1)); dilbert(t); println(t) end
1 Like

Kuramoto model (ten lines)

Kuramoto model is a mathematical model of synchronization. Consider oscillators with random periods. (In tems of the following code, the random periods are denoted by 2\pi/\omega[i] and \omega[i]'s are generated by Normal(1, 2).) The model describes the synchronization of the oscillators:

  • If the strength K of synchronization is less than a certain critical value K_c, then no synchronization will occur.

  • Synchronization begins when the strength K of synchronization exceeds the ctitical value K_c.

The critical value K_c is determined by the probability distribution dist of \omega[i]'s: K_c = π/2/pdf(dist, mean(dist)).

The video below looks like fireflies are synchronizing their rhythmic flashes.

a = 1.3; tmax = 25.0; using Distributions, DifferentialEquations, Plots
kuramoto!(dθ, θ, param, t) = (N = length(θ); (K, ω) = param; 
    for i in 1:N dθ[i] = ω[i] + K*mean(sin(θ[j] - θ[i]) for j in 1:N) end)
(m, n) = (32, 16); (d, v) = (Normal(0, 2), 1.0); K_c = 2/π/pdf(d, 0)
θ₀ = 2π*rand(m, n); tspan = (0.0, tmax); K = a*K_c; ω = rand(d, m, n) .+ v
sol = solve(ODEProblem(kuramoto!, θ₀, tspan, (K, ω)))
anim = @animate for t in [fill(0, 20); 0:0.1:tmax; fill(tmax, 20)]
    plot(size=(400, 220), title="Kuramoto model K = $(a)K_c: t = $t", titlefontsize=10)
    heatmap!(sin.(sol(t))'; c=:bwr, colorbar=false, frame=false, ticks=false)
end; gif(anim, "kuramoto$(a).gif")

Result
kuramoto1.3

For the readability, expand the above compressed code, in an appropriate environment, by

quote
a = 1.3; tmax = 25.0; using Distributions, DifferentialEquations, Plots
kuramoto!(dθ, θ, param, t) = (N = length(θ); (K, ω) = param; 
    for i in 1:N dθ[i] = ω[i] + K*mean(sin(θ[j] - θ[i]) for j in 1:N) end)
(m, n) = (32, 16); (d, v) = (Normal(0, 2), 1.0); K_c = 2/π/pdf(d, 0)
θ₀ = 2π*rand(m, n); tspan = (0.0, tmax); K = a*K_c; ω = rand(d, m, n) .+ v
sol = solve(ODEProblem(kuramoto!, θ₀, tspan, (K, ω)))
anim = @animate for t in [fill(0, 20); 0:0.1:tmax; fill(tmax, 20)]
    plot(size=(400, 220), title="Kuramoto model K = $(a)K_c: t = $t", titlefontsize=10)
    heatmap!(sin.(sol(t))'; c=:bwr, colorbar=false, frame=false, ticks=false)
end; gif(anim, "kuramoto$(a).gif")
end |> Base.remove_linenums! |> x -> display("text/markdown", "```julia\n$x\n```")

Expanded Code

begin
    a = 1.3
    tmax = 25.0
    using Distributions, DifferentialEquations, Plots
    kuramoto!(dθ, θ, param, t) = begin
            N = length(θ)
            (K, ω) = param
            for i = 1:N
                dθ[i] = ω[i] + K * mean((sin(θ[j] - θ[i]) for j = 1:N))
            end
        end
    (m, n) = (32, 16)
    (d, v) = (Normal(0, 2), 1.0)
    K_c = (2 / π) / pdf(d, 0)
    θ₀ = (2π) * rand(m, n)
    tspan = (0.0, tmax)
    K = a * K_c
    ω = rand(d, m, n) .+ v
    sol = solve(ODEProblem(kuramoto!, θ₀, tspan, (K, ω)))
    anim = @animate(for t = [fill(0, 20); 0:0.1:tmax; fill(tmax, 20)]
                plot(size = (400, 220), title = "Kuramoto model K = $(a)K_c: t = $(t)", titlefontsize = 10)
                heatmap!((sin.(sol(t)))'; c = :bwr, colorbar = false, frame = false, ticks = false)
            end)
    gif(anim, "kuramoto$(a).gif")
end
8 Likes

Here’s an electrostatic particle-in-cell code for collisionless plasmas using Fourier fields and nearest grid point particles that squeezes into 7 lines. It’s set up to simulate the two-stream instability:

using FFTW, Plots;N=128;P=64N;dt=1/4N;NT=1024;w=200/P;Δ=1/N;
x=rand(P);v=rand([-1.0,1.0],P);u()=(x.=mod.(x.+v/2*dt,1));
ik=im.*2π*vcat(1,1:N/2,-N/2+1:-1)./N; ic(x)=Int.(ceil.(x));n=zeros(N);
@gif for t in 1:NT; 
  u(); E=real.(ifft((E=fft((n.*=0;for j∈x;n[ic(j*N)]+=w;end;n))./ik;E[1]*=0;E)));
  u(); v+=E[ic(mod.(x*N,N))]*dt;scatter(x,v,ylims=(-3,3));
end every 8

The counter-streaming particles are plotted with velocity against position. Small perturbations to their charge density are unstable. An electrostatic wave, supported by the particles’ charge density, grows exponentially and eventually traps the particles!

7LinePICFourier

P.S. The energy in the electric field and particles are sum(E.^2)*Δ/2 and sum(v.^2/2)*w) respectively.

Edits: collisionless, remove extra ic

30 Likes

One of my so far unachieved goals has been to program an algorithm for PI on my HP calculator. Reading this post reminded me of that goal so I found an algorithm for computing PI and converted it to Julia in 7 lines without a package.

# calculate n iterations of archimedes PI recurrence
function pi_archimedes(n)
    polygon_edge_length_squared = BigFloat(2.0, precision=24)
    for i in 1:n
        polygon_edge_length_squared = 2 - 2 * sqrt(1 - polygon_edge_length_squared / 4)
    end
    return 4 * 2^n * sqrt(polygon_edge_length_squared) / 2
end

To test try this loop.

julia> for n = 0:16
           result = pi_archimedes(n)
           error = result - π
           println("$n iterations $result error $error")
       end

I think it would be straightforward to put the for loop in one line but I find this more readable.

I should be able to program this into my calculator!

3 Likes

I think I found something which could serve if someone wants to know about the problem: https://www.particleincell.com/2015/two-stream-instability/

4 Likes

I can make a self-referential answer :smirk:

using WordCloud
using HTTP
url = "https://discourse.julialang.org/t/seven-lines-of-julia-examples-sought/50416"
wc = url |> HTTP.get |> HTTP.body |> String |> html2text |> wordcloud |> generate!
paint(wc, "seven-lines-of-julia-examples-sought.svg")

22 Likes

This works like a breeze!

Just tried it on the JuliaCon 2021 page JuliaCon 2021 (times are UTC) :: pretalx :

1 Like
using Test
x = 1
y = 2
@test x ≈ y
...
Test Failed
  Expression: x ≈ y
   Evaluated: 1 ≈ 2

Writing tests in julia is actually fun.

  • There is one and only test framework, that everybody uses. In all other languages, I tried there are multiple frameworks in use.
  • It is zero boilerplate and super intuitive. I think I never needed to google how it works and I think in all other languages test frameworks I needed to consult docs/google.
  • It ships with julia, no need to install it.
  • Gives good error messages containing runtime values.
17 Likes

Cool compression of the PIC code. I think that you can remove one of the two ic(x) definitions to make the code even shorter.

Good call :+1: I’ll edit my post. I’ve started a repo called ParticleInCellCodeGolf.jl full of tiny pic codes; only single species, non-relativistic, 1D, electrostatic so far! I’m not sure I can face doing anything more with such a terrible coding style. Everyone should look at GEMPIC.jl to see a proper PIC code.

The world’s simplest climate model. It is not really modelling climate, just plugging in the absorptivity of the sun’s energy coming to earth and the emissivity of the longwave radiation leaving the earth. It does calculate the temperature of the earth as a function of the day of year. It assumes that the earth is isothermal and does not have any thermal mass.
In the line count I exclude the definition of constants. For the sake of readability I have gone two over the number of lines.
And if you don’t like the temperature it is a simple matter to adjust the sun’s output, or the emissivity or absorptivity of the earth! :slight_smile: No need to move!
Another interesting experiment is to equate the emissivity and absorptivity and then adjust them.
With the temperature calculation inside a 4th root the temperature stays fairly stable. (1.7% variation in temperature for a 6.8% variation in heat flux over the course of a year)

#= Calculate the temperature of the earth using the simplest model=#
using Unitful, Plots
p_sun = 386e24u"W"      # power output of the sun
radius_a = 6378u"km"    # semi-major axis of the earth
radius_b = 6357u"km"    # semi-minor axis of the earth
orbit_a = 149.6e6u"km"  # distance from the sun to earth
orbit_e = 0.017         # eccentricity of r = a(1-e^2)/(1+ecos(θ))  & time ≈ 365.25 * θ / 360 where θ is in degrees
a = 0.75                # absorptivity of the sun's radiation
e = 0.6                 # emmissivity of the earth (very dependent on cloud cover)
σ = 5.6703e-8u"W*m^-2*K^-4"  # Stefan-Boltzman constant
temp_sky = 3u"K"        # sky temperature
# constants defined above and not included in the code
t = (0:0.25:365.25)u"d"       # day of year
θ = 2*π/365.25u"d" .* t       # approximate angle around the sun
r = orbit_a * (1-orbit_e^2) ./ (1 .+ orbit_e .* cos.(θ))    # distance from sun to earth
area_projected = π * radius_a * radius_b    # area of earth facing the sun
ec = sqrt(1-radius_b^2/radius_a^2)  # eccentricity of earth
area_surface = 2*π*radius_a^2*(1 + radius_b^2/(ec*radius_b^2)*atanh(ec)) # surface area of the earth
q_in = p_sun * a * area_projected ./ (4 * π .* r.^2) # total heat impacting the earth
temp_earth = (q_in ./ (e*σ*area_surface) .+ temp_sky^4).^0.25 # temperature of the earth
plot(t*u"d^-1", temp_earth*u"K^-1" .- 273.15, label = false, title = "Temperature of Earth", xlabel = "Day", ylabel = "Temperature [C]")

earth_temp

10 Likes

@Jake Cool code (pun intended). Great example of the power of Julia

This function composition trick is handy in economic simulations where x is a vector of prices updated every time period according to a certain rule that may be known to converge.

m = 6
prices = rand(m)
quantities = 1 .+ rand(m)
demand(prices) = 1 ./ prices
update(prices) = prices + 0.5*(demand(prices) .- quantities)
Fⁿ(x, n) = ∘(fill(update, n)...)(x)
@show equilibrium = Fⁿ(prices, 50)
# equilibrium = Fⁿ(prices, 50) = [0.5086155857976737, 0.6478537179305888, 0.7804529267779775, 0.6460790380486849, 0.9611891081266686, 0.6886777705264078]

But is the splatting worth the cost? And is there a way to have it return the intermediate values of x?

I normally use

function Fn(x, n)
    res = Vector{Float64}[x]
    for i in 1:n
        push!(res, update(res[end]))
    end
    return res
end

Another neat use of function composition: Minimizing a convex function f(x) over an affine set Ax = b using a projected Newton’s method and autodiff:

using LinearAlgebra; using ForwardDiff: gradient, hessian                           # Imports
n, m = 30, 10                                                                       # Problem size
A = rand(m, n); b = rand(m)                                                         # Constraint Ax = b
Q = Symmetric(rand(n, n))^2; p = rand(n)                                            # Objective function
f(x) = exp(x' * Q * x + p' * x)
step(x) = x + ([hessian(f, x) A'; A zeros(m, m)]\[-gradient(f, x); zeros(m)])[1:n]  # Newton update step
x_star = ∘(fill(step, 20)...)(A\b)                                                  # Get feas. point, iterate

Proof of optimality: Calculate the dual variables to show that \nabla f(x^*) \in \mathcal{C}(A').

μ = A' \ -gradient(f, x_star)
all(isapprox.(A' * μ + gradient(f, x_star), 0, atol=1e-5))
# true
1 Like

If performance matters you can use foldl as proposed upthread by @FedericoStra:

julia> F(y) = 2y;

julia> Fⁿ(x, n) = foldl((y,_) -> F(y), 1:n; init=x);

julia> Fⁿ(1, 10)
1024

You can use accumulate to get an array with the intermediate values:

julia> Fⁿ′(x, n) = accumulate((y,_) -> F(y), 1:n; init=x);

julia> Fⁿ′(1, 10)
10-element Vector{Int64}:
    2
    4
    8
   16
   32
   64
  128
  256
  512
 1024
7 Likes