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")
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
47 Likes
Plots.jl
Is’nt that from Plots.jl documentation ? It’s still beautifull 
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
Many thanks.
Best wishes,
Niall.