If you're following along, you'll need to checkout the dev branch of Plots.jl. Don't do this if you're running on JuliaBox, you apparently have to run package management commands through the console.
In [1]:
Pkg.checkout("Plots","dev")
In [4]:
using Plots
In [5]:
gr()
Out[5]:
Plots.jl is awesome!Do not use Julia just because you want code to go faster.
Some things will be fast right off the bat, but a lot of that has to do more with solid integration with libraries (linear algebra and FFTs, for instance) than with the compiler innovations.
The developers have spent a lot of time making Julia code go quickly, but well-optimized Julia code, as we'll see in the example below, takes some effort to write.
Random fields show up in a lot of different scientific applications (kriging in geostatistics, for instance), and one might be interested in certain geometric properties of what is known as an excursion set, the set of points for which the field takes values above a certain threshold. Let's see how to simulate some random fields and calculate these geometric properties in Julia.
Gaussian random fields (in which the finite-dimensional distributions are Gaussian) are essentially defined by a covariance function, which describes how the value of the field at two points in space is related. Isotropic and stationary random fields have a covariance function which varies only as a function of distance between the two points, and we'll restrict ourselves to such fields where the covariance function takes the form of a Gaussian.
In [6]:
abstract AbstractCovarianceFunction
In [7]:
abstract StationaryCovarianceFunction <: AbstractCovarianceFunction
In [8]:
Base.call(f::StationaryCovarianceFunction,x,y) = f(sqrt(sumabs2(x-y)))
Out[8]:
In [9]:
doc"""
A Gaussian covariance function is of the form
C(x_1,x_2) = σ^2 e^{-\|x_1-x_2\|^2/2L^2}
where L is the bandwidth of the Gaussian kernel.
"""
type GaussianCovarianceFunction <: StationaryCovarianceFunction
L::Real
σ::Real
end
Base.call(f::GaussianCovarianceFunction,h) = f.σ^2*exp(-h^2/(2*f.L^2))
Out[9]:
In [12]:
?GaussianCovarianceFunction
Out[12]:
In [13]:
gf = GaussianCovarianceFunction(16,1)
plot(x->gf(x),0,64,xlabel="Lag",ylabel="Covariance",legend=false)
Out[13]:
In [14]:
randc(m,n) = randn(m,n)+im*randn(m,n)
function grf(m,n,C::StationaryCovarianceFunction)
tx = collect(0:m-1)
ty = collect(0:n-1)
Rjs = zeros(m,n)
for j in 1:n, i in 1:m
Rjs[i,j] = C([tx[1],ty[1]],[tx[i],ty[j]])
end
Sjs = [Rjs Rjs[:,end] Rjs[:,end:-1:2]]
Sjs2 = [Sjs ; Sjs[end,:] ; Sjs[end:-1:2,1] Sjs[end:-1:2,end:-1:2]]
Γ = fft(Sjs2)
Z = randc(size(Sjs2)...)
fft(sqrt(Γ./(4*m*n)).*Z)[1:m,1:n]
end
Out[14]:
In [15]:
grf(128,128,gf) # Run once for good luck!
@time M = grf(1024,1024,gf);
In [20]:
Z1,Z2 = reim(M)
heatmap(Z1)
Out[20]:
And I said we were interested in the excursion set of a Gaussian random field, so let's see what that looks like:
In [18]:
heatmap(Z1.>1)
Out[18]:
The geometrical properties of excursion sets of random fields are captured by something called the Euler characteristic, $\Gamma$. This essentially measures the number of connected components of the excursion set minus the number of holes in those components. It is useful because we can write down an expectation for it
$$\operatorname{E}[\Gamma(A_u)] = \left(\frac{|A|}{(2\pi)^{3/2}}u +\frac{2|\partial A|}{2\pi}\right)e^{-u^2/2} + \Psi(u)$$
In [23]:
EΓ(x,J) = J^2/(2pi)^(3/2)*x.*exp(-x.^2/2) +
2J/(2pi)*exp(-x.^2/2) +
1/sqrt(2pi)*(1-0.5*(1+erf(x/sqrt(2))))
Out[23]:
And to calculate it, we approximate the excursion set by a lattice of pixels in the image. If two adjacent pixels are both in the excursion set, we connect them with an edge. The approximate Euler characteristic is the number of faces on the resulting graph, minus the number of edges, plus the number of vertices.
In [24]:
function Γ(Z::AbstractMatrix{Float64},u)
m,n = size(Z)
F,E,P = 0,0,0
for j in 1:n-1
for i in 1:m-1
if Z[i,j]>u
P += 1
ns1 = Z[i+1,j]>u
ns2 = Z[i+1,j+1]>u
ns3 = Z[i,j+1]>u
E += ns1+ns2+ns3
F += (ns1&ns2)+(ns2&ns3)
end
end
if Z[m,j]>u
P += 1
E += Z[m,j+1]>u
end
end
for i in 1:m-1
if Z[i,n]>u
P += 1
E += Z[i+1,n]>u
end
end
P += Z[m,n]>u
F-E+P
end
function Γ(Z::AbstractMatrix{Float64},u::AbstractVector)
γs = zeros(Int,length(u))
for i in 1:length(u)
γs[i] = Γ(Z,u[i])
end
γs
end
Out[24]:
In [25]:
plot(x->EΓ(x,64),0,5,legend=false,xlabel="Threshold",ylabel="Gamma")
plot!(x->Γ(Z1,x),0,5)
Out[25]:
The expected Euler characteristic is, after all, a statistic, so let's generate a whole bunch of random fields, calculate the mean Euler characteristic and see how well it matches. And see how quick we can do it.
In [26]:
Zs = Matrix{Float64}[]
@time for i in 1:64
M = grf(512,512,gf)
push!(Zs,real(M))
push!(Zs,imag(M))
end
In [27]:
@time γs = Int[Γ(Z,u) for Z in Zs, u in -5:0.1:5];
In [29]:
γm = vec(mean(γs,1))
γv = vec(std(γs,1));
In [51]:
plot(-5:0.1:5,γm,
xlabel="Threshold",
ylabel="Gamma",
line=:solid,
lab="calculated",
leg=false,
w=3)
plot!(-5:0.1:5,γm+2*sqrt(γv),c=:lightblue,line=:dash,w=3)
plot!(-5:0.1:5,γm-2*sqrt(γv),c=:lightblue,line=:dash,w=3)
plot!(x->EΓ(x,32),-5,5,w=2,lab="theoretical")
Out[51]:
In [ ]: