SPDE menetekel

animate.jl

# ] add https://github.com/mschauer/BridgeSPDE.jl

# will run in the parent directory
# of the script location and create output folder there
cd(joinpath(@__DIR__, ".."))


using SparseArrays
using BridgeSPDE
using FileIO
using Statistics
using LinearAlgebra
using Trajectories
using ImageTransformations
using Random
using Colors
using Makie

inner(x::Vector) = dot(x,x)
inner(x::Vector, y::Vector) = dot(x,y)
inner(x, y) = x'*y
inner(x) = x'*x


println("\ninit")

const F0 = Float32
downscale = 8
ds(img, ds) = imresize(img, floor(Int,size(img,1)/ds), floor(Int,size(img,2)/ds))
pic = "cormullion"
juliapng = FileIO.load(joinpath(@__DIR__, "..", "$pic.png"))
ju  = F0.(2 .- norm.(ds(juliapng, downscale)))[2:end-1,2:end-1]

m1, n1 = size(ju)
pad = 25
julia = zeros(F0, m1 + 2pad, n1 + 2pad)
julia[pad:pad+m1-1, pad:pad+n1-1] = 2ju/3
image(julia)
Random.seed!(1)
θ = F0(0.3), F0(0.3)

mkpath(joinpath(@__DIR__, "..", "output$pic"))

x0 = vec(julia)

m, n = size(julia)

σ = 0.1
dt = F0(0.06)
J = gridderiv(F0, m, n)

J1, J1t = J[1], J[2]
J2, J2t = J[3], J[4]
img(x) = image(reshape(x, (m, n)))
mat(x) = reshape(x, (m, n))
Λ = gridlaplacian(F0, m, n) + 0.01*I

d = m*n

l = 1000
T = sum(l)*dt

θ1, θ2 = θ


t = 0.0
ts = t:dt:T

B = -σ^2/2*Λ + θ1*J1 + θ2*J2t # wind direction/drift north-east
B2 = -σ^2/2*Λ + θ1*J1t + θ2*J2 # adjoint process with wind direction/drift south-west


rs(x) = reshape(x, (m, n))

println("\nsimulate trajectory")
x = copy(x0);
for i in 1:10
    global x
    x = x + dt*(B*x) + σ*sqrt(dt)*randn(F0, d)
end

global X = trajectory((ts[1]=>reshape(x, (m,n)),))
for i in 2:length(ts)
    global x,t
    x = x + dt*(B*x) + σ*sqrt(dt)*randn(F0, d)
    t += dt
    push!(X, ts[i] => reshape(x, (m, n)))

end
reverse!(X.x)
x = X.x[end][:]
for i in 2:length(ts)
    global x,t
    x = x + dt*(B2*x) + σ*sqrt(dt)*randn(F0, d)
    t += dt
    push!(X, ts[i] => reshape(x, (m, n)))
end


@time begin
    println("\ncreate images")

    output_obs = Observable(reshape(clamp.(first(X.x), -1, 1)/2 .+ 0.5, (m, n)))
    scene = surface(output_obs; shading = false, show_axis = false, colormap = :deep)
    scale!(scene, 1.0, 1.0, 17.5)

    st = Stepper(scene, joinpath(dirname(@__DIR__), "output$pic"))

    for (i,j) in enumerate(eachindex(X.x)[1:10:end])
        global output_obs[] = reshape(clamp.(X.x[j], -1, 1)/2 .+ 0.5, (m, n))

        step!(st) # save the current state of the Scene to an image

    end

end

run(`ffmpeg -y -r 40 -f image2 -i output$pic/surf-%d.png -vcodec libx264 -crf 25 -pix_fmt yuv420p surf$pic.mp4`)
#run(`ffmpeg -i surf$pic.mp4 -i music.mp3 -shortest surfm$pic.mp4`)

cormullion.png

see https://chalmersuniversity.box.com/s/j6auchr11k7otpvwii799oi7stgzwvg2

@time begin
    println("\ncreate images")
    
    output_obs = Observable(reshape(clamp.(first(X.x), -1, 1)/2 .+ 0.5, (m, n)))
    scene = surface(output_obs; shading = false, show_axis = false, colormap = :deep)
    scale!(p1, 1.0, 1.0, 17.5)

    st = Stepper(scene, joinpath(dirname(@__DIR__), "output$pic"))

    for (i,j) in enumerate(eachindex(X.x)[1:10:end])
        global output_obs[] = reshape(clamp.(X.x[j], -1, 1)/2 .+ 0.5, (m, n))

        #FileIO.save(joinpath(@__DIR__, "..", "output$pic", "img-$i.png"), Gray.(out))
        
        step!(st) # save the current state of the Scene to an image

    end

end

Thanks!