Seven Lines of Julia (examples sought)

Checked the Laplacian: Δ(u) == ∇⋅∇(u) # and got it is true

5 Likes

Addition of uniform, identically distributed random variables. In one case we are adding a variable to itself (dependent distribution), and in the other we we are adding it to an independent, identically distributed variable.

using Plots, MonteCarloMeasurements

a = 1..2
b = 1..2

default(alpha=0.5, nbins=20, normalize=true)
plot(a+a, label="a+a")
plot!(a+b, label="a+b")

image

8 Likes

Adjacency matrix of hyper-rectangular/toroidal lattices of any dimensions

using LinearAlgebra, SparseArrays
A⊕B = kron(I(size(B,1)), A)+kron(B, I(size(A,1)))
lattice(n) = spdiagm(1=>trues(n-1),-1=>trues(n-1))
lattice(L...) = lattice(L[1])⊕lattice(L[2:end]...)
clattice(n) = (A=lattice(n); A[1,end]=A[end,1]=true; A)
clattice(L...) = clattice(L[1])⊕clattice(L[2:end]...)

I find this pretty damn useful besides being cool

julia> clattice(3,4,2,5) # toroidal lattice of dimensions 3,4,2,5
120×120 SparseMatrixCSC{Int64, Int64} with 840 stored entries:
⢞⣱⢄⠑⠑⢄⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠓⣄⠀⠀⠀⠀⠀⠀
⢄⠑⢞⣱⠀⠀⠑⢄⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠓⣄⠀⠀⠀⠀
⠑⢄⠀⠀⢞⣱⢄⠑⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠓⣄⠀⠀
⠀⠀⠑⢄⢄⠑⢞⣱⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠓⣄
⠑⢄⠀⠀⠀⠀⠀⠀⢞⣱⢄⠑⠑⢄⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠑⢄⠀⠀⠀⠀⢄⠑⢞⣱⠀⠀⠑⢄⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠑⢄⠀⠀⠑⢄⠀⠀⢞⣱⢄⠑⠀⠀⠀⠀⠓⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠑⢄⢄⠑⢞⣱⠀⠀⠀⠀⠀⠀⠓⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⢞⣱⢄⠑⠓⣄⠀⠀⠓⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⢄⠑⢞⣱⠀⠀⠓⣄⠀⠀⠓⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢤⠀⠀⠙⢤⠀⠀⢏⡵⢄⠑⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢤⠀⠀⠙⢤⢄⠑⢏⡵⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢤⠀⠀⠀⠀⠀⠀⢏⡵⢄⠑⠑⢄⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢤⠀⠀⠀⠀⢄⠑⢏⡵⠀⠀⠑⢄⠀⠀⠑⢄⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠑⢄⠀⠀⢏⡵⢄⠑⠀⠀⠀⠀⠑⢄⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠑⢄⢄⠑⢏⡵⠀⠀⠀⠀⠀⠀⠑⢄
⠙⢤⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⢏⡵⢄⠑⠑⢄⠀⠀
⠀⠀⠙⢤⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⢄⠑⢏⡵⠀⠀⠑⢄
⠀⠀⠀⠀⠙⢤⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠑⢄⠀⠀⢏⡵⢄⠑
⠀⠀⠀⠀⠀⠀⠙⢤⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠑⢄⢄⠑⢏⡵
6 Likes

Plotting an animated heart

using Plots
n = 400
t = range(0, 2π, length = n)
x = 16sin.(t).^3
y = 13cos.(t) .- 5cos.(2t) .- 2cos.(3t) .- cos.(4t)
@gif for i ∈ 1:n
    circleplot(x, y, i, line_z = 1:n, cbar = false, c = :reds,framestyle = :none)
end when i > 40 && mod1(i, 10) == 5

jl_K9sme4

47 Likes

Plots.jl

Is’nt that from Plots.jl documentation ? It’s still beautifull :slight_smile:

4 Likes

I get: UndefVarError: circleplot not defined

2 Likes

I think there is no rule the code has to be originally written by the Discourse user. It only need to be short, memorable, and not under a restrictive license.

3 Likes

that is the way of greater benefit

1 Like

My current favourite

using ForwardDiff, ReverseDiff
function callprice(r, T, σ, S₀, K; n = 100)
    memo = Dict(); dt = T / n; u = exp(σ * √dt); q = (exp(r * dt) - 1/u) / (u - 1/u)
    function helper(i, pos)
        get!(memo, (i, pos)) do
            if i == n; max(S₀ * u^pos - K, 0) else exp(- r * dt) * (q*helper(i+1, pos+1) + (1-q)*helper(i+1, pos-1)) end end end
    helper(0, 0) end
(ForwardDiff.derivative(s -> callprice(r, T, σ, s, K), S₀), ReverseDiff.gradient(s -> callprice(r, T, σ, s[1], K), [S₀]))

Automatic differentiation through a memoized recursive function … would probably not even have tried that in any other language.
Unfortunately, I did not get it to work when using an array instead of a dictionary as the memo-table. Maybe I should have another look at Zygote buffers …

Still I would expect mentioning the source - actually this is required by most FOSS licences.

6 Likes

The other day, a friend came to me asking for the survival probabilities of the normal distributions from 9 to 16 sigmas. I answered “Why dont you compute them yourself, you know how they are defined, and you know how to get them in both R and python !”. The precision of the question ticked in my head for a second…

We both know that the Normal distribution is light tailed, so we expected these probabilities to be pretty low. Then, firering up Julia, I understood why he had troubles founding them :

using Distributions
N = Normal();
1 .- cdf(N,1:16)
1 .- cdf(N,big.(1:16))

This is really cool. One of the principal reasons I love Julia.

19 Likes

The survival function/complementary CDF actually comes with Distributions.jl:

using Distributions
N = Normal();
ccdf(N,1:16)
ccdf(N,big.(1:16))
3 Likes

ASCII Mandelbrot set (had to golf it a bit to fit in 7 lines)

function mandelbrot(a, z=0)
    for i=1:2e4
        z = z^2 + a
    end
    abs(z)
end
[prod(x->mandelbrot(complex(x, y)) < 2 ? "∘" : " ",-2.0:0.03125:0.5) for y=-1:0.05:1].|>println;

in the REPL it prints

                                                                ∘                
                                                                                 
                                                                                 
                                                            ∘∘                   
                                                         ∘∘∘∘∘∘∘                 
                                                         ∘∘∘∘∘∘∘                 
                                                          ∘∘∘∘∘                  
                                                            ∘                    
                                               ∘∘    ∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘            
                                               ∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘ ∘∘∘      
                                                ∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘      
                                              ∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘       
                                             ∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘      
                                            ∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘    
                                           ∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘     
                           ∘∘   ∘∘        ∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘    
                            ∘∘∘∘∘∘∘∘∘∘   ∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘    
                         ∘∘∘∘∘∘∘∘∘∘∘∘∘∘  ∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘    
                         ∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘ ∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘     
                     ∘∘∘ ∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘ ∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘      
 ∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘        
                     ∘∘∘ ∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘ ∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘      
                         ∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘ ∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘     
                         ∘∘∘∘∘∘∘∘∘∘∘∘∘∘  ∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘    
                            ∘∘∘∘∘∘∘∘∘∘   ∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘    
                           ∘∘   ∘∘        ∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘    
                                           ∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘     
                                            ∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘    
                                             ∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘      
                                              ∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘       
                                                ∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘      
                                               ∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘ ∘∘∘      
                                               ∘∘    ∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘∘            
                                                            ∘                    
                                                          ∘∘∘∘∘                  
                                                         ∘∘∘∘∘∘∘                 
                                                         ∘∘∘∘∘∘∘                 
                                                            ∘∘                   
                                                                                 
                                                                                 
                                                                ∘ 
28 Likes

I just come across this implementation of Mathematica’s Nest:

nest( f::Function, x0, n::Int) = (n <= 0) ? x0 : (f ∘ nest)(f,x0,n-1)

in my search for an efficient Julia implementation of Mathematica’s NestList. Does anyone have such a thing? Just createing a list recursively like this implementation would be expensive, so it’d be necessary to allocate the list at the beginning then fill it. Is there any way the matrix construction notation [f(x) for i in range] might be useful?

Thanks,
Niall.

1 Like
julia> nestlist(f::Function, x0, n::Integer) = n==0 ? [x0] : push!(nestlist(f, x0, n-1), f(nestlist(f, x0, n-1)[end]))
nestlist (generic function with 1 method)

julia> nestlist(x->x*2, 1, 3)
4-element Vector{Int64}:
 1
 2
 4
 8

This only works if f(x0) has the same type as x0. Otherwise, I need to use an array of elements of generic type, like

nestlist(f::Function, x0, n::Integer) = n==0 ? Any[x0] : push!(nestlist(f, x0, n-1), f(nestlist(f, x0, n-1)[end]))
1 Like

Oh, right - thanks! There are two things that are surprising me here, but I’m relatively new to Julia. First, I had thought that the multiple push!es would bring a significant overhead, and second, I had thought the recursion would have an overhead. So I was thinking more of something along the lines of:

function nestlist( f::Function, x0, n::Integer)
	if n ≤ 0
		[x0]
	else
		nl = fill(x0,n+1)
		for i in 1:n
			nl[i+1] = f(nl[i])
		end
		nl
	end
end

So am I totaly mistaken in these concerns?

1 Like

push! is not expensive as it has amortized cost of O(1), though it’s certainly advantageous to pre-allocate the array and avoid recursion. By the way, if you really want to make it as efficient as possible, you could change nl = fill(x0,n+1) into

nl = Vector{typeof(x0)}(undef, n+1)
nl[1] = x0
2 Likes

The only real performance killer in @greatpet’s code is calling nestlist a second time just to get the last element. And I think this is being done there just to fit it into a single line, instead of two lines.

2 Likes

:grin: Many thanks.

Best wishes,

Niall.