Seven Lines of Julia (examples sought)

For a fairer comparison:

  • Python has mean, median, and geometric_mean functions that could be employed like the Julia library functions.
  • There’s no arbitrary recursion depth here.
5 Likes

Thanks.
It was a snip from the 1st comment I saw on Twitter.

Language syntax comparison suffer from the same type of problem as language performance comparison:
it depends on the ability of the person writing the code & on the actual language

The world would be better off if there was a “language elegance game” just like there is a Benchmarks Game.

7 Likes

:interrobang:Wooaah! remarkable! Could you explain how that magic trick works?

            Fⁿ(x,n)=  ∘(fill(F,n)...)(x) 

Thanks!

1 Like

∘ 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

18 Likes

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)

my_seven_lines

Here’s a bonus animation:
PursuitCurves_Circle_R02

26 Likes

Don’t forget one of the most evocative tools in Julia-land, the splat, which is needed to take the array-typed output of fill and use the elements of that array to populate the individual arguments of .

2 Likes

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.

2021-03-16 19_53_23-Bending-Deflection.pdf - Foxit Reader

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:
my_seven_lines_2_R01

It’s easy enough to change the load function:

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

my_seven_lines_2_R02

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

my_seven_lines_2_R00

17 Likes
  • eight or nine easily read lines are better than seven munched up lines
8 Likes

Hey!
I love your code, but you should use

Fn(n) = reduce(\circ, fill(F, n))(x)
...
Fn(100)[1,1,2,3,5]
1 Like

This is not my code, I saw it on Twitter. But yes it’s true that reducing is faster than splatting. But I think the point here was the clarity of the code not really the performance.

1 Like

Your code still wastes memory to create the useless vector fill(F, n). Also, it doesn’t support the case n=0. Just do this instead

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

Edit: of course one can use a StaticArrays.SVector{3} to sqeeze out the last bit of performance.

3 Likes

Just slightly over, but here’s a user-defined linear regression in Soss.jl:

julia> using Soss, MeasureTheory, SampleChainsDynamicHMC

julia> a = -0.2; b=2.0;

julia> x = randn(100); y = a .+ b .* x;

julia> m = @model x begin
       a ~ Normal()
       b ~ Normal()
       y ~ For(x) do xj Normal(a + b * xj, 1) end
       end;

julia> sample(DynamicHMCChain, m(x=x) | (;y=y))
4000-element MultiChain with 4 chains and schema (b = Float64, a = Float64)
(b = 1.985±0.097, a = -0.207±0.11)
9 Likes

It’s also standard notation for function iteration in the Dynamical Systems literature.
I agree it’s super neat!

Can you help me understand, what does “unsafe” & “type piracy” mean (in general & in this context)?

Btw, for composition I define a Nest(f,x0,n) function like Mathematica:

Nest(f,x0,n) = (n == 0) ? x0 : (f∘Nest)(f,x0,n-1)
import Base.^;
h^(n::Integer) = x0 -> foldl( (x, _) -> h(x), 1:n; init = x0)
2 Likes

Type piracy is when a module defines a method but neither the function nor any of the method’s parameter types are created by that module.

1 Like

Thanks, but what does that mean in this context?

Base.:^(f::Function, n) = n==1 ? f : f ∘ f^(n-1)

That snippet uses ^ and Function, but both of those come from another module, Base. It would not be type piracy to define

struct CustomCounter end

Base.:^(f::Function, n::CustomCounter) = n==1 ? f : f ∘ f^(n-1)

because CustomCounter is defined in this module.

9 Likes

It’s also super slow and wasteful

julia> @time @btime (F^(10^3))([1,1,2,3,5])
  660.172 μs (2001 allocations: 218.89 KiB)
 10.475864 seconds (29.82 M allocations: 3.087 GiB, 4.50% gc time, 0.24% compilation time)
3-element Vector{Float64}:
 2.089057949736859
 2.089057949736859
 2.089057949736859

and doomed to fail miserably

julia> (F^(10^5))([1,1,2,3,5])
...

Bonus: if one is really concerned about the type piracy, then

(f::Function)↑(n::Integer) = x0 -> foldl((x, _) -> f(x), 1:n; init = x0)

is a suggestive alternative ( is \uparrow).

Edit: in the first code snippet, @time @btime is a typo. The intended code was just @time, to take into consideration also the compilation time, instead of @btime.

4 Likes

Try my version w/ ‘foldl’

Yes your version with foldl is perfectly fine!

1 Like

I like Pink noise. It’s is Gaussian noise with power spectral density 1/fᵅ. One can sample it with the inverse Fourier transform


# Make pink noise
using FFTW, Statistics
function pinknoise(m, n) 
    A = real(ifft((Complex{Float64}[i+j==2 ? 0. : (randn()+randn()*im)/sqrt((sin((i-1)*pi/m)^2 + sin((j-1)* pi/n)^2)) for i in 1:m, j in 1:n])))
    real(A)/std(A)
end
A = pinknoise(500, 500)

test

It most spectacular property is that it is invariant under conformal mappings, functions that locally preserves angles, but not necessarily lengths, for example the map x + y*im -> ((x + y*im)/33)^2.

# Make conformal map
getindexfc(A, i, j) = A[mod1(ceil(Int, i), size(A,1)), mod1(ceil(Int, j), size(A,2))]
finv(x) = reim(sqrt(x[1] + x[2]*im)*33) # use inverse for look-up
cfmap(A) = map(CartesianIndices(size(A))) do I
     getindexfc(A, finv(I)...) 
end

This is how it acts on a grid matrix G:

[G  cfmap(G)]

And here with pink noise

[A  cfmap(A)]

PS: I took a bit extra effort to make it tiling, try [A A A; A A A; A A A]

21 Likes