Seven Lines of Julia (examples sought)

Please don’t do it. It has its good and bad sides, but it is fine small snippet, and there are people who like it.

3 Likes

Please please please don’t do it! It doesn’t annoy me, I’m not trying to fight or anything like that. I just wanted to point out the efficiency aspect of it. It can be instructional to understand compile-time/run-time, creation of nested objects with complicated types, etc…

I hope you didn’t take any offence from my posts, and if so I apologize. I wish you a nice weekend :slight_smile:

10 Likes

PS: Anyone interested in camouflage and digital camouflage? I find this interesting from a statistical perspective (you want to match the spectrum of the background) and the mathematical, almost fractal properties of good camouflage.

5 Likes

I used to work with camouflage in cephalopods. It’s a relatively complex field with a lot of interesting opinions… I haven’t published anything recently in the field, and I’m sure a lot of new stuff came along since, but you’d need to catch up on things to make sure you’re not redoing stuff that’s already been published. If you like I can put you in touch with some of the experts in the field.

It’s a cool problem that involves:

  1. detection versus recognition (being “invisible” versus being unrecognizable)
  2. what can the camouflager output in terms of (as you put it) spectrum (how varied is your display) (spatially, chromatically, temporally, etc)
  3. the composition of the background you’ll be viewed on (again, spatially, chromatically, temporally)
  4. what can the visual systems you’re hiding from resolve (yet again, spatially, chromatically, temporally)?

You can see how this becomes very complicated very fast when you consider that everything changes with time of day, day of year, behavior (being motionless, mimicking leaves moving in the wind), and more. The whole distinction between being detected and being recognized makes everything a lot more complicated too, because you can be perfectly detectable (e.g. bright flamboyant colors) but utterly unrecognizable.

7 Likes

I see that here the biologist kicks in and thinks about a catalogue of everything that is out there in nature… I am quite aware of that, but I want to think a bit like a physicist here and impose stronger invariance and symmetry on the problem: e.g. statistical spatial stationarity, scale invariance, or conformal invariance as above. More on Zulip: https://julialang.zulipchat.com/#narrow/stream/225582-random/topic/Camouflage/near/234025423

2 Likes

This isn’t anything out of this world, but it showcases Julia’s ability to represent expressions and to easily perform arbitrary precision arithmetic with GMP/MPFR.

Continued fractions:

function contfrac(x::Real, n::Integer)
    n < 1 && return Int[]
    fpart, ipart = modf(x)
    fpart == 0 ? [Int(ipart)] : [Int(ipart); contfrac(1/fpart, n-1)]
end

foldr((x, y) -> :($x + 1 / $y), contfrac(big(π), 25)) |> println
# 3 + 1 / (7 + 1 / (15 + 1 / (1 + 1 / (292 + 1 / (1 + 1 / (1 + 1 / (1 + 1 / (2 + 1 / (1 + 1 / (3 + 1 / (1 + 1 / (14 + 1 / (3 + 1 / (3 + 1 / (23 + 1 / (1 + 1 / (1 + 1 / (7 + 1 / (4 + 1 / (35 + 1 / (1 + 1 / (1 + 1 / (1 + 1 / 2)))))))))))))))))))))))

foldr((x, y) -> x + 1 // big(y), contfrac(big(π), 75))
# 13926567982265799805873939967043853503//4432964269365850039583559427823398144

If you put the resulting string

foldr((x,y) -> "$x+\\frac1{$y}", contfrac(big(π), 15)) |> println
# 3+\frac1{7+\frac1{15+...}}

in a LaTeX cell you get

3+\frac1{7+\frac1{15+\frac1{1+\frac1{292+\frac1{1+\frac1{1+\frac1{1+\frac1{2+\frac1{1+\frac1{3+\frac1{1+\frac1{14+\frac1{2+\frac1{1}}}}}}}}}}}}}}
Addendum

Of course in general an iterative implementation such as

function contfrac(x::Real, n::Integer)
    cf = Int[]
    for _ in 1:n
        fpart, ipart = modf(x)
        push!(cf, Int(ipart))
        x = 1 / fpart
    end
    cf
end

would be preferred, because it reduces allocations and avoids stack overflows. This was however a bit too long for the seven lines restriction. Notice that in the case of continued fractions, to actually overflow the stack while performing meaningful computations requires working with a quite large precision. The default precision 256 of BigFloat allows to correctly compute only the first 76 terms of the continued fraction of π.

12 Likes

Well you could also use a reduce operation :slight_smile:

contfrac(x, n) = foldl(1:n, init=(Int[], x)) do (l, r), _
       f, i = modf(r)
       return [l; Int(i)], 1/f
end |> first

julia> contfrac(π, 5)
5-element Vector{Int64}:
   3
   7
  15
   1
 292

This could go on one line too (sort of):

contfrac(x, n) = foldl(((l, r), _) -> let (f, i) = modf(r); ([l; Int(i)], 1/f) end, 1:n, init=(Int[], x))[1]

julia> contfrac(π, 5)
5-element Vector{Int64}:
   3
   7
  15
   1
 292
3 Likes

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

22 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
2 Likes

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
15 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

40 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!

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

36 Likes

This works like a breeze!

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

3 Likes
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.
20 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

14 Likes