@show VERSION
@time using DifferentialEquations
@time using Plots
@time using ProgressMeter
@time using Parameters

⪅(x, y) = x < y || x ≈ y

function isin_ngon(n, R, x, y)
    α = 2π/n
    θ = mod(angle(x + im*y), α)
    r = √(x^2 + y^2) / R
    X, Y = r*cos(θ), r*sin(θ)
    (1 - cos(α))*Y ⪅ sin(α)*(1 - X)
end

isin_unit_pentagon(x, y) = isin_ngon(5, 1, x, y)

x = -1:0.1:1
y = -1:0.1:1
D = isin_unit_pentagon.(x, y')
heatmap(x, y, D'; colorbar=false, size=(200, 200))

struct MyParam{TH, TV, TC, TE}
    h⁻²::TH
    V::TV
    C::TC
    E::TE
end

mystep(x::AbstractVector) = x[2] - x[1]
mystep(x::StepRangeLen) = step(x)

function my_param(isin_func, x::AbstractVector, y::AbstractVector)
    @assert mystep(x) == mystep(y)
    h⁻² = 1/mystep(x)^2
    I, J = eachindex(x), eachindex(y)
    D = isin_func.(x, y')
    V = Tuple{Int64, Int64}[]
    C = Tuple{Int64, Int64}[]
    E = [Tuple{Int64, Int64}[] for i in I, j in J]
    for j in J, i in I
        D[i, j] ? push!(V, (i, j)) : push!(C, (i, j))
        for q in (-1, 1), p in (-1, 1)
            i+p in I && j+q in J && D[i+p, j+q] && push!(E[i, j], (i+p, j+q))
        end
    end
    MyParam(h⁻², V, C, E)
end

function fill_complenment!(u, val, p::MyParam)
    @unpack C = p
    for c in C
        u[c...] = val
    end
    u
end

x = -1:0.1:1
y = -1:0.1:1
p = my_param(isin_unit_pentagon, x, y)

@unpack h⁻², V, C, E = p
@show h⁻²
@show V[50]
@show E[V[50]...]
@show C[10];

u0 = zeros(length(x), length(y))
fill_complenment!(u0, NaN, p)'

f!(du, v, u, p, t) = @. du = v

function g!(dv, v, u, p, t)
    @unpack h⁻², V, C, E = p
    @inbounds for v in V
        i, j = v
        es = E[i,j]
        dv[i,j] = h⁻²*(sum(u[e[1], e[2]] for e in es) - length(es)*u[i,j])
    end
end

x = y = range(-1, 1; length=201)
p = my_param(isin_unit_pentagon, x, y)
tspan = (0.0, 10.0)

f(x, y) = exp(-(x^2 + y^2)/2)
u0 = f.(x, y')
v0 = zero(u0)

@time prob = DynamicalODEProblem(g!, f!, v0, u0, tspan, p)
@time sol = solve(prob)
@time sol = solve(prob)
@time sol = solve(prob);

umin = minimum(minimum(u.x[2][v...] for v in p.V) for u in sol.u)
umax = maximum(maximum(u.x[2][v...] for v in p.V) for u in sol.u)

L = 201
ts = range(tspan...; length=L) |> r -> [fill(r[1], 20); r; fill(r[end], 20)]
prog = Progress(length(ts), 0)
@gif for t in ts
    u = fill_complenment!(sol(t).x[2], NaN, p)
    surface(x, y, u'; zlim=(umin, umax), colorbar=false, c=:CMRmap)
    title!("t = $(round(t, digits=1))"; titlefontsize=12)
    next!(prog)
end