Seven Lines of Julia (examples sought)

[jcook]

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!

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

[Jake]

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!

[mschauer]

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

[guoyongzhi]

I can make a self-referential answer

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

[findmyway]

This works like a breeze!

Just tried it on the JuliaCon 2021 page https://pretalx.com/juliacon2021/featured/ :

[jw3126]

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.

[rouckas]

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.

[jcook]

Good call 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.

[Jake]

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! 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]")

[johnh]

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

[maxkapur]

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

[maxkapur]

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

[sijo]

@maxkapur

But is the splatting worth the cost?

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

@maxkapur

And is there a way to have it return the intermediate values of x ?

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

[tdunning]

Here are 6 lines of code (excluding the import and using statements) that define linear, quadratic and cubic Bezier splines using higher order functions.

using Plots
import Base.+, Base.*, Base./

+(f, g) = x -> f(x) + g(x)
*(t::Number, g) = x -> t * g(x)

interpolate(a, b) = t -> (1.0-t)*a + t*b
b1(p1, p2) = interpolate(p1, p2)
b2(p1, p2, p3) = interpolate(b1(p1, p2), b1(p2, p3))
b3(p1, p2, p3, p4) = interpolate(b2(p1, p2, p3), b2(p2, p3, p4))

The idea here is that a Bezier spline of order n is defined by n+1 reference points as an interpolation between the two splines of order n-1 that are defined respectively with the first n points or the last n points. The code renders this definition literally and thus is easy to match up to the intuitive definition.

More importantly, the machine-level code produced by using this very abstract definition is just as compact and efficient as the code produced by an explicit polynomial implementation of a cubic spline. In fact, you can even derive that polynomial definition by giving symbolic values to the spline instead of numerical values. If you do that, the result is not a numerical position, but is the polynomial expressing the spline. You can mix this up by providing numerical values for the points defining the spline, but giving a symbolic value for the interpolating variable t.

As a usage note, to compute the position along a cubic spline defined by points p1, p2, p3, and p4, use b3(p1, p2, p3, p4)(t)(t)(t) where t ranges from 0 to 1. The first application gives you an interpolated quadratic spline. When you feed t to that, you get an interpolated line. When you feed t to that line, you get a point which is your answer. You can use this multiple application to produce pretty visualizations that explain how these splines are constructed, but you mostly would probably rather hide the repeated application in a convenience function.

[RoyiAvital]

This is a great thread.

Might worth contributing to One Liner Hub.
It has no Julia code, so it is time to have some.

[rafael.guerra]

I propose herein an interesting enhancement of the classic moving (rolling or running) average filter (MA filter).

The key thing about the MA filter is that it can be implemented with a simple algorithm that is extremely fast.

That algorithm is written below but with improved boundary conditions. It seems to run faster than the excellent CartesianIndices code given here.

I have discovered that imposing endpoint reflection symmetry does a good job at the extremities. If you know of a reference for this technique, please let me know.

function moving_average(x, w)   # w is MA window half-length
    n, nw = length(x), 2*w + 1
    y = similar(x)
    x0 = @views [reverse(2*x[1] .- x[2:w+1]); x; reverse(2*x[n] .- x[n-w:n-1])]  # endpoint reflection symmetry
    y[1] = c = sum(view(x0,1:nw))
    for i in 2:n; y[i] = c = c + x0[i+nw-1] - x0[i-1]; end
    y/nw
end

One example of the two moving average techniques mentioned:

# m is a odd integer > 1
function moving_average_CI(A::AbstractArray, m::Int)
    out = similar(A)
    R = CartesianIndices(A)
    Ifirst, Ilast = first(R), last(R)
    I1 = m÷2 * oneunit(Ifirst)
    for I in R
        n, s = 0, zero(eltype(out))
        for J in max(Ifirst, I-I1):min(Ilast, I+I1)
            s += A[J]
            n += 1
        end
        out[I] = s/n
    end
    return out
end

using Plots; gr(dpi=600)
t = 0:0.5e-3:0.13;  w = 10;  # time axis and MA window half-length
x = 100*t + 4*sin.(2π*15*t) + 2(rand(length(t)) .- 0.5)
scatter(t, x, c=:blues, ms=2, msw=0.2, label="Input series with random noise", legend=:outertop)
plot!(t, moving_average(x, 10), c=:red, lw=2.2, label="MA with reflection symmetry ($(2*w+1) samples)")
plot!(t, moving_average_CI(x,2*w+1), c=:blue, lw=0.7, label="MA with CartesianIndices ($(2*w+1) samples)")

[goerch]

I’d like to mention @aplavin `s comment Y.A.t.Q : Yet Another @threads Question - #13 by aplavin which looks like the start of a new search engine.

[Bardo]

@rafael.guerra

If you know of a reference for this technique, please let me know.

All reflection references I found were in a density estimation context.
Here you use it as an extrapolation (prediction) method.
Visually appealing, but justified?
BTW, in x0, the endpoints of x are counted twice.

[rafael.guerra]

@Bardo

BTW, in x0, the endpoints of x are counted twice.

Well spotted! x0 has now been improved in code above. And thanks for the references.

[fgerick]

Not seven lines, but I still enjoy it every time I use it:

using Symbolics
@syms x y z
import LinearAlgebra: dot, (×)
const ∂ = Symbolics.derivative
struct ∇; end
∇(ψ::T) where T = [∂(ψ,x), ∂(ψ,y), ∂(ψ,z)]
dot(::Type{∇}, u) = ∂(u[1],x) + ∂(u[2],y) + ∂(u[3],z)
(×)(::Type{∇}, u) = [∂(u[3],y) - ∂(u[2],z), ∂(u[1],z) - ∂(u[3],x), ∂(u[2],x) - ∂(u[1],y)] 
Δ(ψ) = ∂(∂(ψ,x),x) + ∂(∂(ψ,y),y) + ∂(∂(ψ,z),z)
dot(u, ::Type{∇}) = v->map(vi->u[1]*∂(vi,x) + u[2]*∂(vi,y) + u[3]*∂(vi,z), v)

julia> u = [y^2+z*y,x^2+2y,z*x]
3-element Vector{SymbolicUtils.Symbolic{Number}}:
 y^2 + y*z
 x^2 + 2y
 x*z

julia> (u⋅∇)(u)
3-element Vector{Num}:
  (z + 2y)*(x^2 + 2y) + x*y*z
 2(x^2) + 4y + 2x*(y^2 + y*z)
      z*(x^2) + z*(y^2 + y*z)

julia> ∇⋅u
2 + x

julia> ∇×u
3-element Vector{Num}:
           0
       y - z
 2x - z - 2y