Seven Lines of Julia (examples sought)

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)

anim

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.

48 Likes

I love watching particles dance! Awesome Leandro!

1 Like

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

time_animation

21 Likes

@FedericoStra, found interesting to perform a zoom-in loop and recompute your fractal and to use animated gif technique shown by @lmiq to create:
mandelbrot_zooms

xc, yc = -0.55, 0.61; 
x0, x1 = xc - 2, xc + 2; y0, y1 = yc - 2, yc + 2;
anim = @animate for t in 1:50
    x, y = range(x0, x1; length=1000), range(y0, y1; length=1000)
    heatmap(x, y, -log.(mandelbrot.(x' .+ y .* im));aspect_ratio=1,border=:none,legend=:none);
    x0, x1 = (15x0 + x1)/16, (15x1 + x0)/16;  y0, y1 = (15y0 + y1)/16, (15y1 + y0)/16 
end
gif(anim,"mandelbrot_zooms.gif",fps=10)
14 Likes

for reference, what is the extended code that gives your example (a) parallelized (b) with the actual L-J potential (c) both

1 Like

Nice! I was thinking of doing something similar, but I couldn鈥檛 find the time, nor fit it in less than 7 lines (you cheated! :stuck_out_tongue:).

I think that, to make it better, the number of iterations in mandelbrot has to increase the more one zooms in, otherwise the boundary of the fractal becomes too imprecise. Here is a complete version that does so:

using Plots
exprange(start, stop, len) = exp.(range(log(start), log(stop), length=len))
function mandelbrot(z; lim=75) w = z
    for n = 1:lim;  abs2(w) < 4 ? w = w^2 + z : return n end
    lim + 1
end
x鈧, y鈧 = -0.5626805, 0.6422555
anim = @animate for (r, l) in zip(exprange(2, 1.35e-6, 120), exprange(100, 2500, 120))
    x = range(x鈧-r, x鈧+r; length=600); y = range(y鈧-r, y鈧+r; length=600);
    heatmap(x, y, -log.(log.(mandelbrot.(x' .+ y .* im; lim=round(l))));
        legend=:none, border=:none, ticks=:none, size=(600,600), ratio=1)
end
g = gif(anim; fps=12)

Apparently I cannot attach the animated GIF file (which you can find here), so I include only the last frame:

10 Likes

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

Filtered image: filtered

14 Likes

What about a simulation with balls bouncing against each other?
Or a simulation of moving people getting infected with a virus when they meet other people.

I have one of those, but not in a compact code.

From here: https://github.com/m3g/CKP

2 Likes

@LaurentPlagne posted this beautiful example:

See Plot a circle with a given radius with Plots.jl. I don鈥檛 know if he has published the code for it though.

10 Likes
] activate .
instantiate

This is some of my favorite Julia code.

16 Likes

Thank you for the comment :wink:

Unfortunately I don鈥檛 think that it can fit the 7 lines example list (>1000 SLOC)鈥

We use this example to show how to structure a not-too-small project (PkgSkeleton, types, modules, sub-modules, tests, export, dependencies, multiple dispatch, profile, //isme 鈥) in our Julia training lectures.

Although I enjoy reading the examples in this thread, I would second @Tamas_Papp saying that Julia鈥檚 power is more strongly demonstrated through slightly larger examples.

We must take some time to refactor the package from its training material form and publish it : now I have an extra motivation to see if a Julia wizard could compact it significantly :smiley:

6 Likes

Side note: Love the color scheme. What are you using for it?

1 Like

I think seven lines of code is excellent. It showcases what Julia can do, and hopefully shows how Jukia can express scientific notation or equations in a compact way.

However one small pleading from me. Let鈥檚 now go down the road of Obfuscated Perl. Back in the day when Perl was a dominant language for Web programming, and many other things, there were a lot of smart people involved. Still are probably! There were regular 鈥極bfuscated Perl鈥 competitions.
Also the whitespace as code stuff (where you translate code into 鈥rr鈥 meaningful whitespace ) etc.
I know I am being like the Grinch but IMHO these did the language no good

Don鈥檛 get me wrong - there is nothing bad about using a language feature but please, please explain what it does so other people can learn.

3 Likes

I recently started using SatelliteToolbox.jl, so here is my submission. This code calculates the location (latitude, longitude, altitude) of the ISS 鈥渞ight now鈥. The readability suffers a bit trying to stick to seven, but it鈥檚 not too bad. This is based off of the example here (see the link for full details).

using Dates, SatelliteToolbox  
# Download data for ISS
download("https://celestrak.com/NORAD/elements/stations.txt", joinpath(pwd(), "space_stations.txt")) 
# Initialize the orbit propagator
orbp = init_orbit_propagator(Val(:sgp4), read_tle("space_stations.txt")[1])
# Get the current time (Julian date)
rightnow = DatetoJD(now())
# Propogate the orbit to "right now"
o, r_teme, v_teme = propagate_to_epoch!(orbp, rightnow)
# Get the position (radians, radians, meters)
lat, lon, h = ECEFtoGeodetic(rECItoECEF(TEME(), ITRF(), rightnow, get_iers_eop())*r_teme)
# Nice print out after conversions
println("Current location of the ISS: $(rad2deg(lat))掳  $(rad2deg(lon))掳  $(h/1000) km")
6 Likes

@mihalybaci Is Santa tracking data available on 24 December?

:laughing:

Not sure, I poked around here a bit, and didn鈥檛 see anything for Santa. But this site gets its data from NORAD, and they鈥檙e the ones that do the Santa tracking, so it could be hidden in a file somewhere.

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))
鈭侻鈭倁(u,v) = [FDd(u->M(u,v)[1],u), FDd(u->M(u,v)[2],u), FDd(u->M(u,v)[3],u)]
鈭侻鈭倂(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(鈭侻鈭倁(ui, 0.0), 鈭侻鈭倂(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

Mobius_strip_normal

29 Likes

Instead of compressing it further - I鈥檇 be more interested to see it broken out to be more readable :). Awesome plot though!

5 Likes

For what it is worth, Makie鈥檚 藡band藡 function actually makes 2d bands in 3d space鈥 comes in handy for your M枚bius strips:

t = 0:0.1:2pi+0.1
lower = [Point3f0(.2cos(3s),0.0, 0.0) .+ (1 + 0.2sin(3s))*Point3f0(0, cos(s), sin(s)) for s in t]
upper = [Point3f0(-.2cos(3s),0.0, 0.0) .+ (1 - 0.2sin(3s)).*Point3f0(0, cos(s), sin(s)) for s in t]
band(lower, upper, color = [1.0:length(t)梅2; length(t)梅2:-1:1.0; length(t)梅2:-1:1.0; 1.0:length(t)梅2])

12 Likes