text:
### A Pluto.jl notebook ###
# v0.14.8

using Markdown
using InteractiveUtils

# ╔═╡ c14d5f69-8ed3-4ad7-b094-8de94394ad8e
using Revise

# ╔═╡ 3bbf5ec7-c772-4a1a-9eab-c0f26a99e94f
using PolaronMobility

# ╔═╡ f088c85f-4d12-4f4c-8e3e-318652b68a43
using Plots

# ╔═╡ 2b5f21a3-1130-47cd-9938-f49a60048e27
using DataFrames

# ╔═╡ b38c2c58-689b-48f9-b5b1-1499651e1928
using CSV

# ╔═╡ e9685260-9ece-4c3b-b110-096a7c39857b
using QuadGK

# ╔═╡ e495c578-ac0e-46c3-8f2f-aa0541cf3574
using Optim

# ╔═╡ 7f124d90-0430-11ec-129e-bf0ec7dbda06
# ESA Results

# ╔═╡ bf2bb3d0-6f5d-412a-8f9e-b10cb30560bd
begin
	ħ = 1.05457162825e-34
	eV = 1.602176487e-19
	me = 9.10938188e-31
	kB =  1.3806504e-23
	ϵ_0 = 8.854E-12
	amu = 1.660_539_066_60e-27
	Ha = 4.35974820e-18
	Bohr = 5.29177249e-11
	c = 2.99792458e8
end

# ╔═╡ 802b2d3b-d1d8-40c7-9f34-407780385457
α_range = 1:12

# ╔═╡ 30fe032c-f911-4f7e-aa24-1fb473ecf9e0
β_range = vcat([i for i in 0.1:0.1:8.0], [100.0])

# ╔═╡ 71355321-bb21-42d5-9966-60a9438b762c
Ω_range = vcat([1e-200], [i for i in 0.01:0.1:20.01])

# ╔═╡ a1161327-a8e8-44a0-8af7-b14d437ae14c
ω = 1

# ╔═╡ f40b9810-6a0e-4627-a290-80a446d101c8
m = 1

# ╔═╡ 95538cd0-3588-423a-ab1a-95915ee60828
begin
	# Equation 31: The <|X(t) - X(s)|^{-1}> * exp(-|t-w|) effective action.
	A_integrand(x, v, w) = (w^2 * x + (v^2 - w^2) / v * (abs(1 - exp(-v * x))))^(-0.5) * exp(-x)

	A(v, w, α) = π^(-0.5) * α * v * QuadGK.quadgk(x -> A_integrand(x, v, w), BigFloat(0.0), Inf)[1]

	# Equation 33: Lowest Free energy E = -B - A where B = -3/(4v)*(v-w)^2.
	free_energy(v, w, α) = (3 / (4 * v)) * (v - w)^2 - A(v, w, α)
end

# ╔═╡ 5eaec90c-f922-496e-a388-f40a49f38b26
begin
	# Equation 62d in Hellwarth.
	Y(x, v, β) = 1 / (1 - exp(-v * β)) * (1 + exp(-v * β) - exp(-v * x) - exp(v * (x -     β)))

	# Integrand of Equation 62c in Hellwarth.
	A_integrand(x, v, w, β) =  (exp(β - x) + exp(x)) / sqrt(1e-10 + abs(w^2 * x * (1 - x / β)     + Y(x, v, β) * (v^2 - w^2) / v))

	# Equation 62c in Hellwarth.
	A(v, w, α, β) = α * v / (sqrt(π) * (exp(BigFloat(β)) - 1)) * QuadGK.quadgk(x ->       A_integrand(x, v, w, β), BigFloat(0), BigFloat(β / 2))[1]

	# Equation 62b in Hellwarth. Equation 20 in Osaka.
	B(v, w, β) = 3 / β * (log(v / w) - 1 / 2 * log(2 * π * BigFloat(β)) - log(sinh(v *     BigFloat(β) / 2) / sinh(w * BigFloat(β) / 2)))

	# Equation 62e in Hellwarth. Equation 17 in Osaka.
	C(v, w, β) = 3 / 4 * (v^2 - w^2) / v * (coth(v * BigFloat(β) / 2) - 2 / (v * β))

	# Equation 62a in Hellwarth. In paragraph below Equation 22 in Osaka; has extra 1/	  β due to different definition of A, B & C.
	function free_energy(v, w, α, β)
		setprecision(BigFloat, 32)
		a = A(v, w, α, β)
		b = B(v, w, β)
		c = C(v, w, β)
		-(a + b + c)
	end
end

# ╔═╡ 5568b92e-5a20-4771-b440-4de309b6173c
function variation(α, β; v = 0.0, w = 0.0)

    initial = [v, w]

    # Limits of the optimisation.
    lower = [1, 1]
    upper = [200, 200]

    # Osaka Free Energy function to minimise.
    f(x) = free_energy(x[1], x[2], α, β)

    # Use Optim to optimise the free energy function w.r.t v and w.
    solution = Optim.optimize(
        Optim.OnceDifferentiable(f, initial; autodiff = :forward),
        lower,
        upper,
        initial,
        Fminbox(BFGS()),
    )

    # Get v and w values that minimise the free energy.
    v, w = Optim.minimizer(solution)

	println("v: $v, w: $w")

    # Return variational parameters that minimise the free energy.
    return v, w
end

# ╔═╡ 6fc1e642-77ef-49bb-895d-a4006b563f5e
function feynamn_mass(α, v, w, m, ω)
	j(τ) = 1 + (v * ω / w^2) * (1 - w^2 / v^2) * (1 - exp(-v * τ / ω)) / τ
	integrand(τ) = exp(-τ) * τ^(1/2) * j(τ)^(-3/2)
	return m * (1 + α * (v / w)^3 / (3 * π^(1/2)) * quadgk(τ -> integrand(τ), 0.0, Inf)[1])
end

# ╔═╡ ba42efa0-8f26-4ff1-a0ac-a56c626c7525
begin
	v = Array{Float64}(undef, length(β_range), length(α_range))
	w = Array{Float64}(undef, length(β_range), length(α_range))
	FE = Array{Float64}(undef, length(β_range), length(α_range))
	M = Array{Float64}(undef, length(β_range), length(α_range))
	for α in 1:length(α_range), β in 1:length(β_range)
		println("α: $(α_range[α]), β: $(β_range[β])")
		if β == 1
			v[β, α], w[β, α] = variation(α_range[α], β_range[β]; v = 4.0, w = 2.0)
		else
			v[β, α], w[β, α] = variation(α_range[α], β_range[β]; v = v[β-1, α], w = w[β-1, α])
		end
		FE[β, α] = free_energy(v[β, α], w[β, α], α_range[α], β_range[β])
		M[β, α] = feynamn_mass(α_range[α], v[β, α], w[β, α], m, ω)
		CSV.write("v_data.csv", DataFrames.DataFrame([[0.0, [i for i in α_range]...]'; [β_range v]], :auto))
		CSV.write("w_data.csv", DataFrames.DataFrame([[0.0, [i for i in α_range]...]'; [β_range w]], :auto))
		CSV.write("F_data.csv", DataFrames.DataFrame([[0.0, [i for i in α_range]...]'; [β_range FE]], :auto))
		CSV.write("M_data.csv", DataFrames.DataFrame([[0.0, [i for i in α_range]...]'; [β_range M]], :auto))
	end
end

# ╔═╡ c6a84269-a382-4fbd-a102-8a949fe531ec
v

# ╔═╡ 0e9a3b4a-4626-4efb-8925-f493f0ab5f34
w

# ╔═╡ 441cd179-c28f-4e23-ae23-e6a65d7edb46
FE

# ╔═╡ 3b93da85-80a2-466e-b28e-dcf22b4ff0fe
M

# ╔═╡ 2405eb4e-b81c-46fb-94a3-50456b4d7512
begin
	# Integrand of (31) in Feynman I (Feynman 1955, Physical Review, "Slow 	electrons...")
	fF(τ, v, w) = (abs(w^2 * τ + (v^2 - w^2) / v * (1 - exp(- v * τ))))^-0.5 * exp(-τ)

	# (31) in Feynman I
	AF(v,w,α) = π^(-0.5) * α * v * quadgk(τ -> fF(τ, v, w), 0, Inf)[1]

	# (33) in Feynman I
	F(v,w,α) = (3 / (4 * v)) * (v - w)^2 - AF(v, w, α)
end

# ╔═╡ ae352c06-69fe-41fa-ad4c-bafb23235012
function complex_conductivity(Ω, β, α, v, w)

	println("α = $α, Ω = $Ω")

    # FHIP1962, page 1011, eqn (47c).
    R = (v^2 - w^2) / (w^2 * v)

    # FHIP1962, page 1009, eqn (35c).
    D(x) = w^2 / v^2 * (R * (1 - cos(v * x)) * coth(β * v / 2) + x^2 / β - 1im * (R * sin(v * x) + x))

    # FHIP1962, page 1009, eqn (36).
    S(x) = 2 * α / (3 * √π) * (exp(1im * x) + 2 * cos(x) / (exp(β) - 1)) / (D(x))^(3 / 2)

    # FHIP1962, page 1009, eqn (35a).
    integrand(x) = (1 - exp(1im * Ω * x)) * imag(S(x)) / Ω
    1 / (-1im * Ω + 1im * QuadGK.quadgk(x -> integrand(x), 0.0, Inf)[1])
end

# ╔═╡ 6d015de5-3c44-400c-ba60-80d645a2789c
begin
	σ = Array{ComplexF64}(undef, length(Ω_range), length(β_range))
	for α in 1:length(α_range)
		for β in 1:length(β_range), Ω in 1:length(Ω_range)
			println("α: $(α_range[α]), β: $(β_range[β]), Ω: $(Ω_range[Ω])")
			σ[Ω, β] = complex_conductivity(Ω_range[Ω], β_range[β], α_range[α], v[β, α], w[β, α])
		end
		CSV.write("conductivity_data_$α.csv", DataFrames.DataFrame([[0.0, [i for i in β_range]...]'; [Ω_range σ]], :auto))
	end
end

# ╔═╡ 8f052fb6-65bd-4c5f-b931-af616548fb79
σ

# ╔═╡ 38f4dfbc-9a01-494c-8674-0d1b1cc324d3
contour_real = Plots.contourf(Ω_range, β_range, log.(abs.(real.(σ)))', xlabel = "Ω (THz)", ylabel = "T (K)", fill = cgrad(:thermal, rev = false, categorical = true), linewidth = 1.3, color = :grey10, size = (550, 500), tickfontsize = 12, labelfontsize = 12)

# ╔═╡ 22c03bea-a0d5-4291-a568-f29c21d46443
begin
	plot(β_range[1:end-1], v[1:end-1])
	plot!(β_range[1:end-1], w[1:end-1])
end

# ╔═╡ 35d0a3c8-025c-4e63-a65d-0010baf5b61c
plot(β_range[1:end-1], M[1:end-1])

# ╔═╡ Cell order:
# ╠═7f124d90-0430-11ec-129e-bf0ec7dbda06
# ╠═c14d5f69-8ed3-4ad7-b094-8de94394ad8e
# ╠═3bbf5ec7-c772-4a1a-9eab-c0f26a99e94f
# ╠═f088c85f-4d12-4f4c-8e3e-318652b68a43
# ╠═2b5f21a3-1130-47cd-9938-f49a60048e27
# ╠═b38c2c58-689b-48f9-b5b1-1499651e1928
# ╠═e9685260-9ece-4c3b-b110-096a7c39857b
# ╠═e495c578-ac0e-46c3-8f2f-aa0541cf3574
# ╠═bf2bb3d0-6f5d-412a-8f9e-b10cb30560bd
# ╠═802b2d3b-d1d8-40c7-9f34-407780385457
# ╠═30fe032c-f911-4f7e-aa24-1fb473ecf9e0
# ╠═71355321-bb21-42d5-9966-60a9438b762c
# ╠═a1161327-a8e8-44a0-8af7-b14d437ae14c
# ╠═f40b9810-6a0e-4627-a290-80a446d101c8
# ╠═95538cd0-3588-423a-ab1a-95915ee60828
# ╠═5568b92e-5a20-4771-b440-4de309b6173c
# ╠═5eaec90c-f922-496e-a388-f40a49f38b26
# ╠═6fc1e642-77ef-49bb-895d-a4006b563f5e
# ╠═ba42efa0-8f26-4ff1-a0ac-a56c626c7525
# ╠═c6a84269-a382-4fbd-a102-8a949fe531ec
# ╠═0e9a3b4a-4626-4efb-8925-f493f0ab5f34
# ╠═441cd179-c28f-4e23-ae23-e6a65d7edb46
# ╠═3b93da85-80a2-466e-b28e-dcf22b4ff0fe
# ╠═2405eb4e-b81c-46fb-94a3-50456b4d7512
# ╠═ae352c06-69fe-41fa-ad4c-bafb23235012
# ╠═6d015de5-3c44-400c-ba60-80d645a2789c
# ╠═8f052fb6-65bd-4c5f-b931-af616548fb79
# ╠═38f4dfbc-9a01-494c-8674-0d1b1cc324d3
# ╠═22c03bea-a0d5-4291-a568-f29c21d46443
# ╠═35d0a3c8-025c-4e63-a65d-0010baf5b61c


text:
using ApproxFun, MultivariateOrthogonalPolynomials, Plots, SparseArrays
import ApproxFun: Vec, PiecewiseSegment, ZeroOperator, Block
    import MultivariateOrthogonalPolynomials: DirichletTriangle

d = [Triangle(Vec(1,1),Vec(0,1),Vec(0,2)) , Triangle(Vec(1,2),Vec(0,2),Vec(1,1)),
    Triangle(Vec(0,0), Vec(1,0), Vec(0,1)) , Triangle(Vec(1,1),Vec(1,0),Vec(0,1)) ,
    Triangle(Vec(1,0), Vec(2,0), Vec(1,1)) , Triangle(Vec(2,1), Vec(1,1), Vec(2,0))]

p = plot()
    for ▴ in d
        plot!(▴)
    end
    p

∂d = components(PiecewiseSegment([Vec(0,2), Vec(0,1), Vec(0,0), Vec(1,0), Vec(2,0), Vec(2,1), Vec(1,1), Vec(1,2), Vec(0,2)]))

for s in ∂d
    plot!(s)
end
p

# ιd = [Segment(Vec(0,2), Vec(1,1)), Segment(Vec(0,1), Vec(1,1)),
#         Segment(Vec(0,1),Vec(1,0)), Segment(Vec(1,0), Vec(1,1)),
#         Segment(Vec(1,1), Vec(2,0))]

# interfaces
ιd = Dict{NTuple{2,Int}, Segment{Vec{2,Int}}}()
    ιd[(1,2)] = Segment(Vec(0,2), Vec(1,1))
    ιd[(1,4)] = Segment(Vec(0,1), Vec(1,1))
    ιd[(3,4)] = Segment(Vec(0,1),Vec(1,0))
    ιd[(4,5)] = Segment(Vec(1,0), Vec(1,1))
    ιd[(5,6)] = Segment(Vec(1,1), Vec(2,0))

keys(ιd)

ds = vcat(fill.(DirichletTriangle{1,1,1}.(d),3)...) # repeat each triangle 3 times
rs = [Legendre.(∂d)...,
        Legendre(ιd[1,4]), Legendre(ιd[4,5]),  # straight interface
        Legendre(ιd[1,2]), Legendre(ιd[3,4]), Legendre(ιd[5,6]),
        Legendre(ιd[1,4]), Legendre(ιd[4,5]),  # straight interface
        Legendre(ιd[1,2]), Legendre(ιd[3,4]), Legendre(ιd[5,6]),
        vcat(fill.(JacobiTriangle.(d),3)...)...] # diagonal interface

# rs = [Legendre.(∂d); fill.(Legendre.(values(ιd)),2)...; fill.(JacobiTriangle.(d),3)...]

ui  = T -> 1 + (T-1)*3
ux = T -> 2 + (T-1)*3
uy = T -> 3 + (T-1)*3

Dx = Derivative([1,0])
Dy = Derivative([0,1])

N = 30
sprs = A -> (global N; sparse(A[Block.(1:N), Block.(1:N)]))

A = Matrix{SparseMatrixCSC{Float64,Int}}(undef, length(rs), length(ds))
    for K=1:length(rs), J=1:length(ds) # fill with zeros
        A[K,J] = ZeroOperator(ds[J],rs[K]) |> sprs
    end
    # add boundary conditions
    K = 0;
    K +=1; T = 1; J = ui(T); A[K,J] = (I : ds[J] → rs[K]) |> sprs
    K +=1; T = 3; J = ui(T); A[K,J] = (I : ds[J] → rs[K]) |> sprs
    K +=1; T = 3; J = ui(T); A[K,J] = (I : ds[J] → rs[K]) |> sprs
    K +=1; T = 5; J = ui(T); A[K,J] = (I : ds[J] → rs[K]) |> sprs
    K +=1; T = 6; J = ui(T); A[K,J] = (I : ds[J] → rs[K]) |> sprs
    K +=1; T = 6; J = ui(T); A[K,J] = (I : ds[J] → rs[K]) |> sprs
    K +=1; T = 2; J = ui(T); A[K,J] = (I : ds[J] → rs[K]) |> sprs
    K +=1; T = 2; J = ui(T); A[K,J] = (I : ds[J] → rs[K]) |> sprs
    # add dirichlet interface conditions
    for (T1,T2) in ((1,4), (4,5), (1,2), (3,4), (5,6))
        global K +=1;
        J = ui(T1); A[K,J] = (I : ds[J] → rs[K]) |> sprs
        J = ui(T2); A[K,J] = -(I : ds[J] → rs[K]) |> sprs
    end
    # add lr neumann
    K +=1;
    T = 1; J = uy(T); A[K,J] = (I : ds[J] → rs[K]) |> sprs
    T = 4; J = uy(T); A[K,J] = -(I : ds[J] → rs[K]) |> sprs
    # add ud neumann
    K +=1;
    T = 4; J = ux(T); A[K,J] = (I : ds[J] → rs[K]) |> sprs
    T = 5; J = ux(T); A[K,J] = -(I : ds[J] → rs[K]) |> sprs
    # add diagonal neumann
    for (T1,T2) in ((1,2), (3,4), (5,6))
        global K +=1;
        J = ux(T1); A[K,J] = (I : ds[J] → rs[K]) |> sprs
        J = uy(T1); A[K,J] = (I : ds[J] → rs[K]) |> sprs
        J = ux(T2); A[K,J] = -(I : ds[J] → rs[K]) |> sprs
        J = uy(T2); A[K,J] = -(I : ds[J] → rs[K]) |> sprs
    end

    for T in 1:length(d)
        @show K,T
        global K +=1;
        J = ui(T); A[K,J] = (Dx : ds[J] → rs[K]) |> sprs
        J = ux(T); A[K,J] = -(I : ds[J] → rs[K]) |> sprs
        global K +=1;
        J = ui(T); A[K,J] = (Dy : ds[J] → rs[K]) |> sprs
        J = uy(T); A[K,J] = -(I : ds[J] → rs[K]) |> sprs
        global K +=1;
        J = ux(T); A[K,J] = (Dx : ds[J] → rs[K]) |> sprs
        J = uy(T); A[K,J] = (Dy : ds[J] → rs[K]) |> sprs
    end

M = hvcat(ntuple(_ -> size(A,2),size(A,1)), permutedims(A)...)

rhs = vcat(coefficients.(Fun.(Ref((x,y) -> real(exp(x+im*y))), rs[1:length(∂d)], N))...)

rhs = vcat(pad.(coefficients.(Fun.(Ref((x,y) -> x^2), rs[1:length(∂d)])),N)...)

F = factorize(M)
u_cfs = F \ pad(rhs, size(M,1))

u1 = Fun(ds[4], u_cfs[(4-1)*sum(1:N)+1:4*sum(1:N)])
u1(0.99,0.99)
plot(abs.(u_cfs); yscale=:log10)

σ = svdvals(Matrix(M))
plot(svdvals(Matrix(M)))

u1.coefficients
u1(0.1,1.2)-real(exp(0.1+im*1.2))

plot(abs.([norm((M*u_cfs - pad(rhs, size(M,1)))[N*(K-1)+1:N*K]) for K=1:length(rs)] ); yscale=:log10)

length(rs)
K = 30; norm((M*u_cfs - pad(rhs, size(M,1)))[N*(K-1)+1:N*K])

U = Fun.(Ref((x,y) -> real(exp(x+im*y))), d, sum(1:N))

u_cfs = Vector{Float64}()
    for T in d
        append!(u_cfs, pad(Fun((x,y) -> exp(x)*cos(y), DirichletTriangle{1,1,1}(T)).coefficients, sum(1:N)))
        append!(u_cfs, pad(Fun((x,y) -> exp(x)*cos(y), DirichletTriangle{1,1,1}(T)).coefficients, sum(1:N)))
        append!(u_cfs, pad(Fun((x,y) -> -exp(x)*sin(y), DirichletTriangle{1,1,1}(T)).coefficients, sum(1:N)))
    end

u_cfs

((M*u_cfs) - pad(rhs, size(M,1))) |> norm
((M*u_cfs) - pad(rhs, size(M,1)))[1:(N*(length(rs)-6*3))] |> norm
NN = (N*(length(rs)-6*3))
((M*u_cfs) - pad(rhs, size(M,1)))[NN+1:NN+sum(1:N)] |> norm
((M*u_cfs) - pad(rhs, size(M,1)))[NN+sum(1:N):NN+2sum(1:N)] |> norm

T = d[1]
u1 = pad(Fun((x,y) -> exp(x)*cos(y), DirichletTriangle{1,1,1}(T)).coefficients, sum(1:N))
u1y = pad(Fun((x,y) -> -exp(x)*sin(y), DirichletTriangle{1,1,1}(T)).coefficients, sum(1:N))

A[19,1]*u1 + A[19,3]*u1y

NN

size(M)
size(u_cfs)
size(rhs)

M*u_cfs - pad(rhs,size(M,1))  |> norm

(Dy*u1)(0.1,1.2)

Fun((x,y) -> -exp(x)*sin(y), DirichletTriangle{1,1,1}(T))(0.1,1.2)

A[18,1]*u_cfs[1:sum(1:N)] +
    A[18,2]*u_cfs[sum(1:N)+1:2sum(1:N)]
å
(length(rs)-6*3)

length(rs)

6*

(I : ds[1] → rs[1])[Block.(1:N), Block.(1:N)] * u_cfs[1:210]

M[1:20,1:210]u_cfs[1:210]

Matrix(M[1:20,211:end]) |> norm

A[1,:]

u1.coefficients

## Interlace operator

M

N,M = length(rs), length(ds)
    A = Matrix{Operator{Float64}}(undef, N,M)
    for K=1:N, J=1:M # fill with zeros
        A[K,J] = ZeroOperator(ds[J],rs[K])
    end
    # add boundary conditions
    K = 0;
    K +=1; T = 1; J = ui(T); A[K,J] = I : ds[J] → rs[K]
    K +=1; T = 3; J = ui(T); A[K,J] = I : ds[J] → rs[K]
    K +=1; T = 3; J = ui(T); A[K,J] = I : ds[J] → rs[K]
    K +=1; T = 5; J = ui(T); A[K,J] = I : ds[J] → rs[K]
    K +=1; T = 6; J = ui(T); A[K,J] = I : ds[J] → rs[K]
    K +=1; T = 6; J = ui(T); A[K,J] = I : ds[J] → rs[K]
    K +=1; T = 2; J = ui(T); A[K,J] = I : ds[J] → rs[K]
    K +=1; T = 2; J = ui(T); A[K,J] = I : ds[J] → rs[K]
    # add dirichlet interface conditions
    for (T1,T2) in ((1,4), (4,5), (1,2), (3,4), (5,6))
        global K +=1;
        J = ui(T1); A[K,J] = I : ds[J] → rs[K]
        J = ui(T2); A[K,J] = -I : ds[J] → rs[K]
    end
    # add lr neumann
    K +=1;
    T = 1; J = uy(T); A[K,J] = I : ds[J] → rs[K]
    T = 4; J = uy(T); A[K,J] = -I : ds[J] → rs[K]
    # add up neumann
    K +=1;
    T = 4; J = ux(T); A[K,J] = I : ds[J] → rs[K]
    T = 5; J = ux(T); A[K,J] = -I : ds[J] → rs[K]
    # add diagonal neumann
    for (T1,T2) in ((1,2), (3,4), (5,6))
        global K +=1;
        J = ux(T1); A[K,J] = I : ds[J] → rs[K]
        J = uy(T1); A[K,J] = I : ds[J] → rs[K]
        J = ux(T2); A[K,J] = -I : ds[J] → rs[K]
        J = uy(T2); A[K,J] = -I : ds[J] → rs[K]
    end

    for T in 1:length(d)
        global K +=1;
        J = ui(T); A[K,J] = I : ds[J] → rs[K]
        J = ux(T); A[K,J] = -Dx : ds[J] → rs[K]
        global K +=1;
        J = ui(T); A[K,J] = I : ds[J] → rs[K]
        J = uy(T); A[K,J] = -Dy : ds[J] → rs[K]
        global K +=1;
        J = ux(T); A[K,J] = Dx : ds[J] → rs[K]
        J = uy(T); A[K,J] = Dy : ds[J] → rs[K]
    end

L = Operator(A)

N = 20
M = sparse(L[Block.(1:N), Block.(1:N)])

rhs = Fun.(Ref((x,y) -> real(exp(x+im*y))), rs[1:length(∂d)], N)

u_cfs = M \ pad(vcat(cfs...), size(M,1))

u_cfs[1:sum(1:N)]

Fun(rhs)

rangespace(L)

J
ds[7]

[randn(2,2) randn(2,2);
 randn(2,2) randn(2,2)]

U = Fun.(Ref((x,y) -> real(exp(x+im*y))), d, sum(1:N))
    u_cfs = Vector{Float64}()
    for u in U
        append!(u_cfs, u.coefficients)
        append!(u_cfs, (Dx*u).coefficients)
        append!(u_cfs, (Dy*u).coefficients)
    end

M*u_cfs
ds[1]
rs[1]

(I : ds[1] → rs[1])*U[1]

# add boundary conditions
for K = 1:length(∂d)
    for J = 1:length(d)
        if domain(rs[K]) ⊆ domain(ds[J])
            A[K,J] = I : ds[J] → rs[K]
            break
        end
    end
end

Ñ = length(∂d)
for (K, J1, J2) in ((Ñ+1, 1, 4),
                    (Ñ+2, 1, 2))
    A[K,J1] =  I : ds[J1] → rs[K]
    A[K,J2] = -I : ds[J2] → rs[K]
end

Operator(A)[1:20,1:20]

rs

# add continuity
for K = 1:length(ιd)

∂d[1] ⊆ d[1]

rs[3]

typeof(components(d))

length(ιd)

Vec{2,Float64}
Segment(Vec(0,0), Vec(1.0,0)) |> typeof

import Makie

S = DirichletTriangle{1,1,1}.(components(d))

C1 = I : S[1] → Legendre(Segment(Vec(0,0),Vec(1,0)))
C2 = I : S[1] → Legendre(Segment(Vec(1,0),Vec(0,0)))

a = Legendre(Segment(Vec(0,0),Vec(1,0)))
b = Legendre(Segment(Vec(1,0),Vec(0,0)))

@which Conversion(a,b)

DirichletTriangle{0,1,1}(d[1])

2
C1.op.op.ops

C2.op.op.ops

C1 == C2

C1[1:20,1:20] == C2[1:20,1:20]

Legendre(Segment(Vec(0,0),Vec(1,0)))

∂d = PiecewiseSegment([Vec(0.,0), Vec(1.,0), Vec(1,1), Vec(1,2), Vec(0,1), Vec(-1,1.5), Vec(0,0)])

o = Fun.(S, Ref([1.0]))

o[1]

ext: jl

max_stars_repo_path: examples/L-shaped Laplace.jl

max_issues_repo_path: examples/L-shaped Laplace.jl

max_forks_repo_path: examples/L-shaped Laplace.jl


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:light
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.4'
#       jupytext_version: 1.2.4
#   kernelspec:
#     display_name: Julia 1.3.0-rc4
#     language: julia
#     name: julia-1.3
# ---

# # Access Argo data using [ArgoVis](https://argovis.colorado.edu/ng/home) in Julia
#
# Citation for the Argovis web application and the Argovis database:
# Tucker, T., D. Giglio, M. Scanderbeg, and S.S.P. Shen, 0: Argovis: A Web Application for Fast Delivery, Visualization, and Analysis of Argo Data. J. Atmos. Oceanic Technol., 0, https://doi.org/10.1175/JTECH-D-19-0041.1
#
# Argo data reference
# " These data were collected and made freely available by the International Argo Program and the national programs that contribute to it. (http://www.argo.ucsd.edu, http://argo.jcommops.org). The Argo Program is part of the Global Ocean Observing System. "
# Argo (2000). Argo float data and metadata from Global Data Assembly Centre (Argo GDAC). SEANOE. http://doi.org/10.17882/42182

# +
#run(pipeline(`which python`,"whichpython.txt")) #external python path
#ENV["PYTHON"]=readline("whichpython.txt")
#import Pkg; Pkg.build("PyCall")

using PyCall

requests = pyimport("requests")
pd = pyimport("pandas")
np = pyimport("numpy")
sp = pyimport("scipy.interpolate")
griddata=sp.griddata
datetime = pyimport("datetime")
pdb = pyimport("pdb")
os = pyimport("os")
netCDF4 = pyimport("netCDF4")
netcdf_dataset=netCDF4.Dataset

# +
function get_monthly_profile_pos(month, year)
    baseURL = "https://argovis.colorado.edu/selection/profiles"
    url = baseURL*"/"*string(month)*"/"*string(year)
    resp = requests.get(url)
    monthlyProfilePos = resp.json()
    return monthlyProfilePos
end

function parse_meta_into_df(profiles)
    #initialize dict
    df = pd.DataFrame(profiles)
    return df
end

# set month and year for metadata
month = 1
year = 2004
metaProfiles = get_monthly_profile_pos(month, year)
metaDf = parse_meta_into_df(metaProfiles)
# -


text:
function _generate_julia_command(use_sysimage::Bool)
    if use_sysimage
        return "/opt/bin/julia -J/opt/simplecontainergenerator_containers/sysimage/SimpleContainerGeneratorSysimage.so"
    else
        return "/opt/bin/julia"
    end
end

function _generate_do_not_use_sysimage_julia_script_content(config::Config = Config())
    return _generate_julia_script_content(false, config)
end

function _generate_use_sysimage_julia_script_content(config::Config = Config())
    make_sysimage = config.make_sysimage
    return _generate_julia_script_content(make_sysimage, config)
end

function _generate_julia_script_content(use_sysimage::Bool,
                                                        config::Config = Config())
    julia_command = _generate_julia_command(use_sysimage)
    wrapper_script_env_vars = config.wrapper_script_env_vars
    set_wrapper_script_env_vars_string = ""
    for (name, value) in wrapper_script_env_vars
        set_wrapper_script_env_vars_string *= "export $(name)=\"$(value)\"\n"
    end
    lines = String[
        "$(set_wrapper_script_env_vars_string)",
        "export JULIA_DEPOT_PATH=\"/opt/simplecontainergenerator_containers/julia_depot\"",
        "export JULIA_PROJECT=\"/opt/simplecontainergenerator_containers/julia_project\"",
        "$(julia_command) \"\$@\"",
    ]
    content = join(lines, "\n") * "\n"
    return content
end

function _generate_global_startup_file_content(config::Config = Config())
    lines = String[
        "_simplecontainergenerator_containers_temp_depot_tmpdir = mktempdir()",
        "atexit(() -> rm(_simplecontainergenerator_containers_temp_depot_tmpdir; force = true, recursive = true))",
        "_simplecontainergenerator_containers_temp_depot = joinpath(_simplecontainergenerator_containers_temp_depot_tmpdir, \"simplecontainergenerator_containers_temp_depot\")",
        "mkpath(_simplecontainergenerator_containers_temp_depot)",
        "pushfirst!(Base.DEPOT_PATH, \"/opt/simplecontainergenerator_containers/julia_depot\")",
        "if get(ENV, \"SIMPLECONTAINERGENERATOR_CONTAINER_NO_TEMP_DEPOT\", \"\") != \"true\"",
        "pushfirst!(Base.DEPOT_PATH, _simplecontainergenerator_containers_temp_depot)",
        "end",
    ]
    content = join(lines, "\n") * "\n"
    return content
end


text:
# included in module Wut

INITIAL_BUCKETS = 1000  # must be an even positive integer

mutable struct RandomDistributedScalarEncoder  <: AbstractEncoder

    w::Int32             # width of the contiguous 1 bits, must be odd
    minIndex::Int32
    maxIndex::Int32
    offset::Float64
    oidx::Int32
    n::Int32             # output bit width (must be greater than w)
    resolution::Float64  # input value difference >= resolution do not give identical bits
    name::AbstractString # optional descriptive string
    verbosity::Int32
    maxoverlap::Int32
    random::MersenneTwister
    bucketmap::Dict{Int32,Vector{Int32}}
    numTries::Int32
    maxBuckets::Int32

    function RandomDistributedScalarEncoder(;resolution=0.0,w=21,n=400,name="",offset=NaN,seed=42,verbosity=0)

        minIndex = maxIndex = INITIAL_BUCKETS/2
        mao=2
        bm = Dict{Int32,Vector{Int32}}()
        rng=MersenneTwister(seed)
        t=shuffle(Vector(1:n))
        bm[minIndex]=t[1:w]
        nm = (name != "") ? name : @sprintf("[%s]" , resolution)

        new(w,minIndex,maxIndex,offset,0,n,resolution,nm,verbosity,mao,rng,bm,0,INITIAL_BUCKETS)
    end

end
export RandomDistributedScalarEncoder

function getWidth(e::RandomDistributedScalarEncoder)
    e.n
end
export getWidth

function getDescription(e::RandomDistributedScalarEncoder)
    e.name
end
export getDescription

function getName(e::RandomDistributedScalarEncoder)
    e.name
end
export getName

function repOK(e::RandomDistributedScalarEncoder,v::Vector)
   true
end

function newRepresentation(e::RandomDistributedScalarEncoder, nidx::Int32, idx::Int32)

    t=shuffle(Vector(1:e.n))
    x=shuffle(e.bucketmap[idx])
    n=vcat(t[1],x[2:e.w])
    while ! repOK(e,n)
        t=shuffle(Vector(1:e.n))
        x=shuffle(e.bucketmap[idx])
        n=vcat(t[1],x[2:e.w])
    end
    n

end

function createBucket(e::RandomDistributedScalarEncoder, idx::Int32)

    z=convert(Int32,0)
    if idx < e.minIndex
        z=e.minIndex
        if idx == e.minIndex - 1
            e.minIndex-=1
        else
            createBucket(e,convert(Int32,idx+1))
        end
    else
        z=e.maxIndex
        if idx == e.maxIndex + 1
            e.maxIndex+=1
        else
            createBucket(e,convert(Int32,idx-1))
        end
    end

    e.bucketmap[idx]=newRepresentation(e,idx,z)

end

function mapBucketIndexToNonZeroBits(e::RandomDistributedScalarEncoder, index::Int32)

    idx = index < 1 ? 1 : ( index > e.maxBuckets ? e.maxBuckets : index )

    if ! haskey(e.bucketmap,idx)
        createBucket(e,idx)
    end

    e.bucketmap[idx]
end

function getBucketIndices(e::RandomDistributedScalarEncoder,x::Float64)

    if isnan(e.offset)
        e.offset=convert(Float64,x)
    end

    bucketIdx = convert(Int32,((e.maxBuckets/2) + convert(Int32,round((x - e.offset) / e.resolution))))

    if bucketIdx < 1
      bucketIdx = 1
    elseif bucketIdx >= e.maxBuckets
      bucketIdx = e.maxBuckets
    end
    bucketIdx
end
export getBucketIndices

function encodeIntoArray(e::RandomDistributedScalarEncoder,n::Number,b::BitPat;learn=true, offset=0, length=0)
    fill!(b.b,false)

    idx=getBucketIndices(e,convert(Float64,n))

    nzb=mapBucketIndexToNonZeroBits(e,idx)

    for i in nzb
        b.b[i]=true
    end

    b
end
export encodeIntoArray

#function decode(e::RandomDistributedScalarEncoder, b::BitPat; parentFieldName="")
#
#end
#export decode

function getBucketValues(e::RandomDistributedScalarEncoder)

end
export getBucketValues

function getBucketInfo(e::RandomDistributedScalarEncoder)

end
export getBucketInfo

function topDownCompute(e::RandomDistributedScalarEncoder)

end
export topDownCompute

function closenessScores(e::RandomDistributedScalarEncoder)

end
export closenessScores

# end # module


text:
module UnaryOps
    using ..SuiteSparseGraphBLAS: isGxB, isGrB, TypedUnaryOperator, AbstractUnaryOp
    using ..libgb
    export UnaryOp
    function UnaryOp(name)
        if isGxB(name) || isGrB(name) #If it's a GrB/GxB op we don't want the prefix
            simplifiedname = name[5:end]
        else
            simplifiedname = name
        end
        tname = Symbol(simplifiedname * "_T")
        simplifiedname = Symbol(simplifiedname)
        structquote = quote
            struct $tname <: AbstractUnaryOp
                typedops::Dict{DataType, TypedUnaryOperator}
                name::String
                $tname() = new(Dict{DataType, TypedUnaryOperator}(), $name)
            end
        end
        @eval($structquote) #Eval the struct into the Types submodule to avoid clutter.
        constquote = quote
            const $simplifiedname = $tname()
            export $simplifiedname
        end
        @eval($constquote)
        return getproperty(UnaryOps, simplifiedname)
    end

    struct GenericUnaryOp <: AbstractUnaryOp
        typedops::Dict{DataType, TypedUnaryOperator}
        name::String
        GenericUnaryOp(name) = new(Dict{DataType, TypedUnaryOperator}(), name)
    end
end

const UnaryUnion = Union{AbstractUnaryOp, TypedUnaryOperator}

#TODO: Rewrite
function _createunaryops()
    builtins = [
    "GrB_IDENTITY",
    "GrB_AINV",
    "GxB_LNOT",
    "GrB_MINV",
    "GxB_ONE",
    "GrB_ABS",
    "GrB_BNOT",
    "GxB_SQRT",
    "GxB_LOG",
    "GxB_EXP",
    "GxB_LOG2",
    "GxB_SIN",
    "GxB_COS",
    "GxB_TAN",
    "GxB_ACOS",
    "GxB_ASIN",
    "GxB_ATAN",
    "GxB_SINH",
    "GxB_COSH",
    "GxB_TANH",
    "GxB_ASINH",
    "GxB_ACOSH",
    "GxB_ATANH",
    "GxB_SIGNUM",
    "GxB_CEIL",
    "GxB_FLOOR",
    "GxB_ROUND",
    "GxB_TRUNC",
    "GxB_EXP2",
    "GxB_EXPM1",
    "GxB_LOG10",
    "GxB_LOG1P",
    "GxB_LGAMMA",
    "GxB_TGAMMA",
    "GxB_ERF",
    "GxB_ERFC",
    "GxB_FREXPE",
    "GxB_FREXPX",
    "GxB_CONJ",
    "GxB_CREAL",
    "GxB_CIMAG",
    "GxB_CARG",
    "GxB_ISINF",
    "GxB_ISNAN",
    "GxB_ISFINITE",
    "GxB_POSITIONI",
    "GxB_POSITIONI1",
    "GxB_POSITIONJ",
    "GxB_POSITIONJ1",
]
    for name ∈ builtins
        UnaryOps.UnaryOp(name)
    end
end

#This is adapted from the fork by cvdlab.
#Add a new GrB_UnaryOp to an AbstractUnaryOp.
function _addunaryop(op::AbstractUnaryOp, fn::Function, ztype::GBType{T}, xtype::GBType{U}) where {T, U}
    function unaryopfn(z, x)
        unsafe_store!(z, fn(x))
        return nothing
    end
    opref = Ref{libgb.GrB_UnaryOp}()
    unaryopfn_C = @cfunction($unaryopfn, Cvoid, (Ptr{T}, Ref{U}))
    libgb.GB_UnaryOp_new(opref, unaryopfn_C, ztype, xtype, op.name)
    op.typedops[U] = TypedUnaryOperator{xtype, ztype}(opref[])
    return nothing
end

#UnaryOp constructors
#####################
function UnaryOps.UnaryOp(name::String, fn::Function, ztype, xtype; keep=false)
    @warn "Use built-in functions where possible, user defined functions are less performant."
    length(name) == 0 && (name = string(fn))
    if keep
        op = UnaryOps.UnaryOp(name)
    else
        op = UnaryOps.GenericUnaryOp(name)
    end
    _addunaryop(op, fn, toGBType(ztype), toGBType(xtype))
    return op
end
#Same xtype, ztype.
function UnaryOps.UnaryOp(name::String, fn::Function, type; keep=false)
    return UnaryOps.UnaryOp(name, fn, type, type; keep)
end
#Vector of xtypes and ztypes, add a GrB_UnaryOp for each.
function UnaryOps.UnaryOp(name::String, fn::Function, ztype::Vector{DataType}, xtype::Vector{DataType}; keep=false)
    @warn "Use built-in functions where possible, user defined functions are less performant."
    length(name) == 0 && (name = string(fn))
    if keep
        op = UnaryOps.UnaryOp(name)
    else
        op = UnaryOps.GenericUnaryOp(name)
    end
    length(ztype) == length(xtype) || throw(DimensionMismatch("Lengths of ztype and xtype must match."))
    for i ∈ 1:length(ztype)
        _addunaryop(op, fn, toGBType(ztype[i]), toGBType(xtype[i]))
    end
    return op
end
#Vector but same ztype xtype.
function UnaryOps.UnaryOp(name::String, fn::Function, type::Vector{DataType}; keep=false)
    return UnaryOps.UnaryOp(name, fn, type, type; keep)
end
#Construct it using all the built in primitives.
function UnaryOps.UnaryOp(name::String, fn::Function; keep=false)
    return UnaryOps.UnaryOp(name, fn, valid_vec; keep)
end

function UnaryOps.UnaryOp(fn::Function; keep=false)
    return UnaryOps.UnaryOp("", fn; keep)
end

function _load(unaryop::AbstractUnaryOp)
    booleans = ["GrB_IDENTITY", "GrB_AINV", "GrB_MINV", "GxB_LNOT", "GxB_ONE", "GrB_ABS"]
    integers = [
        "GrB_IDENTITY",
        "GrB_AINV",
        "GrB_MINV",
        "GxB_LNOT",
        "GxB_ONE",
        "GrB_ABS",
        "GrB_BNOT",
    ]
    unsignedintegers = [
        "GrB_IDENTITY",
        "GrB_AINV",
        "GrB_MINV",
        "GxB_LNOT",
        "GxB_ONE",
        "GrB_ABS",
        "GrB_BNOT",
    ]
    floats = [
        "GrB_IDENTITY",
        "GrB_AINV",
        "GrB_MINV",
        "GxB_LNOT",
        "GxB_ONE",
        "GrB_ABS",
        "GxB_SQRT",
        "GxB_LOG",
        "GxB_EXP",
        "GxB_LOG2",
        "GxB_SIN",
        "GxB_COS",
        "GxB_TAN",
        "GxB_ACOS",
        "GxB_ASIN",
        "GxB_ATAN",
        "GxB_SINH",
        "GxB_COSH",
        "GxB_TANH",
        "GxB_ASINH",
        "GxB_ACOSH",
        "GxB_ATANH",
        "GxB_SIGNUM",
        "GxB_CEIL",
        "GxB_FLOOR",
        "GxB_ROUND",
        "GxB_TRUNC",
        "GxB_EXP2",
        "GxB_EXPM1",
        "GxB_LOG10",
        "GxB_LOG1P",
        "GxB_LGAMMA",
        "GxB_TGAMMA",
        "GxB_ERF",
        "GxB_ERFC",
        "GxB_FREXPE",
        "GxB_FREXPX",
        "GxB_ISINF",
        "GxB_ISNAN",
        "GxB_ISFINITE",
    ]
    positionals = ["GxB_POSITIONI", "GxB_POSITIONI1", "GxB_POSITIONJ", "GxB_POSITIONJ1"]
    complexes = [
        "GxB_IDENTITY",
        "GxB_AINV",
        "GxB_MINV",
        "GxB_ONE",
        "GxB_SQRT",
        "GxB_LOG",
        "GxB_EXP",
        "GxB_LOG2",
        "GxB_SIN",
        "GxB_COS",
        "GxB_TAN",
        "GxB_ACOS",
        "GxB_ASIN",
        "GxB_ATAN",
        "GxB_SINH",
        "GxB_COSH",
        "GxB_TANH",
        "GxB_ASINH",
        "GxB_ACOSH",
        "GxB_ATANH",
        "GxB_SIGNUM",
        "GxB_CEIL",
        "GxB_FLOOR",
        "GxB_ROUND",
        "GxB_TRUNC",
        "GxB_EXP2",
        "GxB_EXPM1",
        "GxB_LOG10",
        "GxB_LOG1P",
        "GxB_CONJ",
        "GxB_CREAL",
        "GxB_CIMAG",
        "GxB_CARG",
        "GxB_ABS",
        "GxB_ISINF",
        "GxB_ISNAN",
        "GxB_ISFINITE",
    ]
    name = unaryop.name
    if name ∈ booleans
        constname = name * "_BOOL"
        unaryop.typedops[Bool] = TypedUnaryOperator(load_global(constname, libgb.GrB_UnaryOp))
    end

    if name ∈ integers
        unaryop.typedops[Int8] = TypedUnaryOperator(load_global(name * "_INT8", libgb.GrB_UnaryOp))
        unaryop.typedops[Int16] = TypedUnaryOperator(load_global(name * "_INT16", libgb.GrB_UnaryOp))
        unaryop.typedops[Int32] = TypedUnaryOperator(load_global(name * "_INT32", libgb.GrB_UnaryOp))
        unaryop.typedops[Int64] = TypedUnaryOperator(load_global(name * "_INT64", libgb.GrB_UnaryOp))
    end

    if name ∈ unsignedintegers
        unaryop.typedops[UInt8] = TypedUnaryOperator(load_global(name * "_UINT8", libgb.GrB_UnaryOp))
        unaryop.typedops[UInt16] = TypedUnaryOperator(load_global(name * "_UINT16", libgb.GrB_UnaryOp))
        unaryop.typedops[UInt32] = TypedUnaryOperator(load_global(name * "_UINT32", libgb.GrB_UnaryOp))
        unaryop.typedops[UInt64] = TypedUnaryOperator(load_global(name * "_UINT64", libgb.GrB_UnaryOp))
    end

    if name ∈ floats
        unaryop.typedops[Float32] = TypedUnaryOperator(load_global(name * "_FP32", libgb.GrB_UnaryOp))
        unaryop.typedops[Float64] = TypedUnaryOperator(load_global(name * "_FP64", libgb.GrB_UnaryOp))
    end
    if name ∈ positionals
        unaryop.typedops[Any] = TypedUnaryOperator(load_global(name * "_INT64", libgb.GrB_UnaryOp))
    end
    name = "GxB_" * name[5:end]
    if name ∈ complexes
        unaryop.typedops[ComplexF32] = TypedUnaryOperator(load_global(name * "_FC32", libgb.GrB_UnaryOp))
        unaryop.typedops[ComplexF64] = TypedUnaryOperator(load_global(name * "_FC64", libgb.GrB_UnaryOp))
    end
end

ztype(::TypedUnaryOperator{I, O}) where {I, O} = O
xtype(::TypedUnaryOperator{I, O}) where {I, O} = I


text:
# Julia wrapper for header: /usr/local/include/sodium.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

@c Int32 sodium_init () libsodium
@c Int32 crypto_auth_hmacsha512256 (Ptr{Uint8}, Ptr{Uint8}, Uint64, Ptr{Uint8}) libsodium
@c Int32 crypto_auth_hmacsha512256_verify (Ptr{Uint8}, Ptr{Uint8}, Uint64, Ptr{Uint8}) libsodium
@c size_t crypto_auth_bytes () libsodium
@c size_t crypto_auth_keybytes () libsodium
@c Ptr{Uint8} crypto_auth_primitive () libsodium
@c Int32 crypto_auth (Ptr{Uint8}, Ptr{Uint8}, Uint64, Ptr{Uint8}) libsodium
@c Int32 crypto_auth_verify (Ptr{Uint8}, Ptr{Uint8}, Uint64, Ptr{Uint8}) libsodium
@c Int32 crypto_auth_hmacsha256 (Ptr{Uint8}, Ptr{Uint8}, Uint64, Ptr{Uint8}) libsodium
@c Int32 crypto_auth_hmacsha256_verify (Ptr{Uint8}, Ptr{Uint8}, Uint64, Ptr{Uint8}) libsodium
@c Int32 crypto_box_curve25519xsalsa20poly1305 (Ptr{Uint8}, Ptr{Uint8}, Uint64, Ptr{Uint8}, Ptr{Uint8}, Ptr{Uint8}) libsodium
@c Int32 crypto_box_curve25519xsalsa20poly1305_open (Ptr{Uint8}, Ptr{Uint8}, Uint64, Ptr{Uint8}, Ptr{Uint8}, Ptr{Uint8}) libsodium
@c Int32 crypto_box_curve25519xsalsa20poly1305_keypair (Ptr{Uint8}, Ptr{Uint8}) libsodium
@c Int32 crypto_box_curve25519xsalsa20poly1305_beforenm (Ptr{Uint8}, Ptr{Uint8}, Ptr{Uint8}) libsodium
@c Int32 crypto_box_curve25519xsalsa20poly1305_afternm (Ptr{Uint8}, Ptr{Uint8}, Uint64, Ptr{Uint8}, Ptr{Uint8}) libsodium
@c Int32 crypto_box_curve25519xsalsa20poly1305_open_afternm (Ptr{Uint8}, Ptr{Uint8}, Uint64, Ptr{Uint8}, Ptr{Uint8}) libsodium
@c size_t crypto_box_publickeybytes () libsodium
@c size_t crypto_box_secretkeybytes () libsodium
@c size_t crypto_box_beforenmbytes () libsodium
@c size_t crypto_box_noncebytes () libsodium
@c size_t crypto_box_zerobytes () libsodium
@c size_t crypto_box_boxzerobytes () libsodium
@c size_t crypto_box_macbytes () libsodium
@c Ptr{Uint8} crypto_box_primitive () libsodium
@c Int32 crypto_box_keypair (Ptr{Uint8}, Ptr{Uint8}) libsodium
@c Int32 crypto_box_beforenm (Ptr{Uint8}, Ptr{Uint8}, Ptr{Uint8}) libsodium
@c Int32 crypto_box_afternm (Ptr{Uint8}, Ptr{Uint8}, Uint64, Ptr{Uint8}, Ptr{Uint8}) libsodium
@c Int32 crypto_box_open_afternm (Ptr{Uint8}, Ptr{Uint8}, Uint64, Ptr{Uint8}, Ptr{Uint8}) libsodium
@c Int32 crypto_box (Ptr{Uint8}, Ptr{Uint8}, Uint64, Ptr{Uint8}, Ptr{Uint8}, Ptr{Uint8}) libsodium
@c Int32 crypto_box_open (Ptr{Uint8}, Ptr{Uint8}, Uint64, Ptr{Uint8}, Ptr{Uint8}, Ptr{Uint8}) libsodium
@c Int32 crypto_core_hsalsa20 (Ptr{Uint8}, Ptr{Uint8}, Ptr{Uint8}, Ptr{Uint8}) libsodium
@c Int32 crypto_core_salsa20 (Ptr{Uint8}, Ptr{Uint8}, Ptr{Uint8}, Ptr{Uint8}) libsodium
@c Int32 crypto_core_salsa2012 (Ptr{Uint8}, Ptr{Uint8}, Ptr{Uint8}, Ptr{Uint8}) libsodium
@c Int32 crypto_core_salsa208 (Ptr{Uint8}, Ptr{Uint8}, Ptr{Uint8}, Ptr{Uint8}) libsodium
@c Int32 crypto_generichash_blake2b (Ptr{Uint8}, size_t, Ptr{Uint8}, Uint64, Ptr{Uint8}, size_t) libsodium
@c Int32 crypto_generichash_blake2b_init (Ptr{crypto_generichash_blake2b_state}, Ptr{Uint8}, size_t, size_t) libsodium
@c Int32 crypto_generichash_blake2b_update (Ptr{crypto_generichash_blake2b_state}, Ptr{Uint8}, Uint64) libsodium
@c Int32 crypto_generichash_blake2b_final (Ptr{crypto_generichash_blake2b_state}, Ptr{Uint8}, size_t) libsodium
@c size_t crypto_generichash_bytes_min () libsodium
@c size_t crypto_generichash_bytes_max () libsodium
@c size_t crypto_generichash_keybytes_min () libsodium
@c size_t crypto_generichash_keybytes_max () libsodium
@c size_t crypto_generichash_blockbytes () libsodium
@c Ptr{Uint8} crypto_generichash_primitive () libsodium
@c Int32 crypto_generichash (Ptr{Uint8}, size_t, Ptr{Uint8}, Uint64, Ptr{Uint8}, size_t) libsodium
@c Int32 crypto_generichash_init (Ptr{crypto_generichash_blake2b_state}, Ptr{Uint8}, size_t, size_t) libsodium
@c Int32 crypto_generichash_update (Ptr{crypto_generichash_blake2b_state}, Ptr{Uint8}, Uint64) libsodium
@c Int32 crypto_generichash_final (Ptr{crypto_generichash_blake2b_state}, Ptr{Uint8}, size_t) libsodium
@c Int32 crypto_hash_sha512 (Ptr{Uint8}, Ptr{Uint8}, Uint64) libsodium
@c Int32 crypto_hash (Ptr{Uint8}, Ptr{Uint8}, Uint64) libsodium
@c Int32 crypto_hash_sha256 (Ptr{Uint8}, Ptr{Uint8}, Uint64) libsodium
@c Int32 crypto_hashblocks_sha256 (Ptr{Uint8}, Ptr{Uint8}, Uint64) libsodium
@c Int32 crypto_hashblocks_sha512 (Ptr{Uint8}, Ptr{Uint8}, Uint64) libsodium
@c Ptr{Uint8} crypto_onetimeauth_poly1305_ref_implementation_name () libsodium
@c Int32 crypto_onetimeauth_poly1305_set_implementation (Ptr{crypto_onetimeauth_poly1305_implementation},) libsodium
@c Ptr{crypto_onetimeauth_poly1305_implementation} crypto_onetimeauth_pick_best_implementation () libsodium
@c Int32 crypto_onetimeauth_poly1305 (Ptr{Uint8}, Ptr{Uint8}, Uint64, Ptr{Uint8}) libsodium
@c Int32 crypto_onetimeauth_poly1305_verify (Ptr{Uint8}, Ptr{Uint8}, Uint64, Ptr{Uint8}) libsodium
@c size_t crypto_onetimeauth_bytes () libsodium
@c size_t crypto_onetimeauth_keybytes () libsodium
@c Ptr{Uint8} crypto_onetimeauth_primitive () libsodium
@c Int32 crypto_onetimeauth (Ptr{Uint8}, Ptr{Uint8}, Uint64, Ptr{Uint8}) libsodium
@c Int32 crypto_onetimeauth_verify (Ptr{Uint8}, Ptr{Uint8}, Uint64, Ptr{Uint8}) libsodium
@c Int32 crypto_scalarmult_curve25519 (Ptr{Uint8}, Ptr{Uint8}, Ptr{Uint8}) libsodium
@c Int32 crypto_scalarmult_curve25519_base (Ptr{Uint8}, Ptr{Uint8}) libsodium
@c size_t crypto_scalarmult_bytes () libsodium
@c size_t crypto_scalarmult_scalarbytes () libsodium
@c Ptr{Uint8} crypto_scalarmult_primitive () libsodium
@c Int32 crypto_scalarmult_base (Ptr{Uint8}, Ptr{Uint8}) libsodium
@c Int32 crypto_scalarmult (Ptr{Uint8}, Ptr{Uint8}, Ptr{Uint8}) libsodium
@c Int32 crypto_secretbox_xsalsa20poly1305 (Ptr{Uint8}, Ptr{Uint8}, Uint64, Ptr{Uint8}, Ptr{Uint8}) libsodium
@c Int32 crypto_secretbox_xsalsa20poly1305_open (Ptr{Uint8}, Ptr{Uint8}, Uint64, Ptr{Uint8}, Ptr{Uint8}) libsodium
@c size_t crypto_secretbox_keybytes () libsodium
@c size_t crypto_secretbox_noncebytes () libsodium
@c size_t crypto_secretbox_zerobytes () libsodium
@c size_t crypto_secretbox_boxzerobytes () libsodium
@c Ptr{Uint8} crypto_secretbox_primitive () libsodium
@c Int32 crypto_secretbox (Ptr{Uint8}, Ptr{Uint8}, Uint64, Ptr{Uint8}, Ptr{Uint8}) libsodium
@c Int32 crypto_secretbox_open (Ptr{Uint8}, Ptr{Uint8}, Uint64, Ptr{Uint8}, Ptr{Uint8}) libsodium
@c Int32 crypto_shorthash_siphash24 (Ptr{Uint8}, Ptr{Uint8}, Uint64, Ptr{Uint8}) libsodium
@c size_t crypto_shorthash_bytes () libsodium
@c size_t crypto_shorthash_keybytes () libsodium
@c Ptr{Uint8} crypto_shorthash_primitive () libsodium
@c Int32 crypto_shorthash (Ptr{Uint8}, Ptr{Uint8}, Uint64, Ptr{Uint8}) libsodium
@c Int32 crypto_sign_ed25519 (Ptr{Uint8}, Ptr{Uint64}, Ptr{Uint8}, Uint64, Ptr{Uint8}) libsodium
@c Int32 crypto_sign_ed25519_open (Ptr{Uint8}, Ptr{Uint64}, Ptr{Uint8}, Uint64, Ptr{Uint8}) libsodium
@c Int32 crypto_sign_ed25519_keypair (Ptr{Uint8}, Ptr{Uint8}) libsodium
@c Int32 crypto_sign_ed25519_seed_keypair (Ptr{Uint8}, Ptr{Uint8}, Ptr{Uint8}) libsodium
@c size_t crypto_sign_bytes () libsodium
@c size_t crypto_sign_publickeybytes () libsodium
@c size_t crypto_sign_secretkeybytes () libsodium
@c Ptr{Uint8} crypto_sign_primitive () libsodium
@c Int32 crypto_sign_seed_keypair (Ptr{Uint8}, Ptr{Uint8}, Ptr{Uint8}) libsodium
@c Int32 crypto_sign_keypair (Ptr{Uint8}, Ptr{Uint8}) libsodium
@c Int32 crypto_sign (Ptr{Uint8}, Ptr{Uint64}, Ptr{Uint8}, Uint64, Ptr{Uint8}) libsodium
@c Int32 crypto_sign_open (Ptr{Uint8}, Ptr{Uint64}, Ptr{Uint8}, Uint64, Ptr{Uint8}) libsodium
@c Int32 crypto_sign_edwards25519sha512batch (Ptr{Uint8}, Ptr{Uint64}, Ptr{Uint8}, Uint64, Ptr{Uint8}) libsodium
@c Int32 crypto_sign_edwards25519sha512batch_open (Ptr{Uint8}, Ptr{Uint64}, Ptr{Uint8}, Uint64, Ptr{Uint8}) libsodium
@c Int32 crypto_sign_edwards25519sha512batch_keypair (Ptr{Uint8}, Ptr{Uint8}) libsodium
@c Int32 crypto_stream_xsalsa20 (Ptr{Uint8}, Uint64, Ptr{Uint8}, Ptr{Uint8}) libsodium
@c Int32 crypto_stream_xsalsa20_xor (Ptr{Uint8}, Ptr{Uint8}, Uint64, Ptr{Uint8}, Ptr{Uint8}) libsodium
@c Int32 crypto_stream_xsalsa20_beforenm (Ptr{Uint8}, Ptr{Uint8}) libsodium
@c Int32 crypto_stream_xsalsa20_afternm (Ptr{Uint8}, Uint64, Ptr{Uint8}, Ptr{Uint8}) libsodium
@c Int32 crypto_stream_xsalsa20_xor_afternm (Ptr{Uint8}, Ptr{Uint8}, Uint64, Ptr{Uint8}, Ptr{Uint8}) libsodium
@c size_t crypto_stream_keybytes () libsodium
@c size_t crypto_stream_noncebytes () libsodium
@c Ptr{Uint8} crypto_stream_primitive () libsodium
@c Int32 crypto_stream (Ptr{Uint8}, Uint64, Ptr{Uint8}, Ptr{Uint8}) libsodium
@c Int32 crypto_stream_xor (Ptr{Uint8}, Ptr{Uint8}, Uint64, Ptr{Uint8}, Ptr{Uint8}) libsodium
@c Int32 crypto_stream_aes128ctr (Ptr{Uint8}, Uint64, Ptr{Uint8}, Ptr{Uint8}) libsodium
@c Int32 crypto_stream_aes128ctr_xor (Ptr{Uint8}, Ptr{Uint8}, Uint64, Ptr{Uint8}, Ptr{Uint8}) libsodium
@c Int32 crypto_stream_aes128ctr_beforenm (Ptr{Uint8}, Ptr{Uint8}) libsodium
@c Int32 crypto_stream_aes128ctr_afternm (Ptr{Uint8}, Uint64, Ptr{Uint8}, Ptr{Uint8}) libsodium
@c Int32 crypto_stream_aes128ctr_xor_afternm (Ptr{Uint8}, Ptr{Uint8}, Uint64, Ptr{Uint8}, Ptr{Uint8}) libsodium
@c Int32 crypto_stream_salsa20 (Ptr{Uint8}, Uint64, Ptr{Uint8}, Ptr{Uint8}) libsodium
@c Int32 crypto_stream_salsa20_xor (Ptr{Uint8}, Ptr{Uint8}, Uint64, Ptr{Uint8}, Ptr{Uint8}) libsodium
@c Int32 crypto_stream_salsa20_beforenm (Ptr{Uint8}, Ptr{Uint8}) libsodium
@c Int32 crypto_stream_salsa20_afternm (Ptr{Uint8}, Uint64, Ptr{Uint8}, Ptr{Uint8}) libsodium
@c Int32 crypto_stream_salsa20_xor_afternm (Ptr{Uint8}, Ptr{Uint8}, Uint64, Ptr{Uint8}, Ptr{Uint8}) libsodium
@c Int32 crypto_stream_salsa2012 (Ptr{Uint8}, Uint64, Ptr{Uint8}, Ptr{Uint8}) libsodium
@c Int32 crypto_stream_salsa2012_xor (Ptr{Uint8}, Ptr{Uint8}, Uint64, Ptr{Uint8}, Ptr{Uint8}) libsodium
@c Int32 crypto_stream_salsa2012_beforenm (Ptr{Uint8}, Ptr{Uint8}) libsodium
@c Int32 crypto_stream_salsa2012_afternm (Ptr{Uint8}, Uint64, Ptr{Uint8}, Ptr{Uint8}) libsodium
@c Int32 crypto_stream_salsa2012_xor_afternm (Ptr{Uint8}, Ptr{Uint8}, Uint64, Ptr{Uint8}, Ptr{Uint8}) libsodium
@c Int32 crypto_stream_salsa208 (Ptr{Uint8}, Uint64, Ptr{Uint8}, Ptr{Uint8}) libsodium
@c Int32 crypto_stream_salsa208_xor (Ptr{Uint8}, Ptr{Uint8}, Uint64, Ptr{Uint8}, Ptr{Uint8}) libsodium
@c Int32 crypto_stream_salsa208_beforenm (Ptr{Uint8}, Ptr{Uint8}) libsodium
@c Int32 crypto_stream_salsa208_afternm (Ptr{Uint8}, Uint64, Ptr{Uint8}, Ptr{Uint8}) libsodium
@c Int32 crypto_stream_salsa208_xor_afternm (Ptr{Uint8}, Ptr{Uint8}, Uint64, Ptr{Uint8}, Ptr{Uint8}) libsodium
@c Int32 crypto_verify_16 (Ptr{Uint8}, Ptr{Uint8}) libsodium
@c Int32 crypto_verify_32 (Ptr{Uint8}, Ptr{Uint8}) libsodium
@c Int32 randombytes_set_implementation (Ptr{randombytes_implementation},) libsodium
@c None randombytes (Ptr{Uint8}, Uint64) libsodium
@c Ptr{Uint8} randombytes_implementation_name () libsodium
@c uint32_t randombytes_random () libsodium
@c None randombytes_stir () libsodium
@c uint32_t randombytes_uniform (uint32_t,) libsodium
@c None randombytes_buf (Ptr{Void}, size_t) libsodium
@c Int32 randombytes_close () libsodium
@c Ptr{Uint8} randombytes_salsa20_implementation_name () libsodium
@c uint32_t randombytes_salsa20_random () libsodium
@c None randombytes_salsa20_random_stir () libsodium
@c uint32_t randombytes_salsa20_random_uniform (uint32_t,) libsodium
@c None randombytes_salsa20_random_buf (Ptr{Void}, size_t) libsodium
@c Int32 randombytes_salsa20_random_close () libsodium
@c Ptr{Uint8} randombytes_sysrandom_implementation_name () libsodium
@c uint32_t randombytes_sysrandom () libsodium
@c None randombytes_sysrandom_stir () libsodium
@c uint32_t randombytes_sysrandom_uniform (uint32_t,) libsodium
@c None randombytes_sysrandom_buf (Ptr{Void}, size_t) libsodium
@c Int32 randombytes_sysrandom_close () libsodium
@c Ptr{Uint8} _sodium_alignedcalloc (Ptr{Ptr{Uint8}}, size_t) libsodium
@c None sodium_memzero (Ptr{Void}, size_t) libsodium
@c Int32 sodium_memcmp (Ptr{Void}, Ptr{Void}, size_t) libsodium
@c Ptr{Uint8} sodium_version_string () libsodium
@c Int32 sodium_library_version_major () libsodium
@c Int32 sodium_library_version_minor () libsodium


text:
"""
    kmeans(X::AbstractMatrix{<:Real}, r::AbstractVector{Int},
        c::AbstractVector{Int}, μ::AbstractMatrix{<:Real}; <keyword arguments>)

Cluster the data `X` with the ``k``-means algorithm.

# Keyword arguments
- `maxiter::Int=256`: the number of maximum iterations.
- `dist::SemiMetric=SqEuclidean()`: the distance function.
- `ϵ::AbstractFloat=1.0e-6`: the absolute tolerance for convergence.
"""
function kmeans(X::AbstractMatrix{<:Real}, r::AbstractVector{Int},
                c::AbstractVector{Int}, μ::AbstractMatrix{<:Real};
                maxiter::Int=256, dist::SemiMetric=SqEuclidean(),
                ϵ::AbstractFloat=1.0e-6)
    mx, n = size(X)
    mμ, k = size(μ)

    mx == mμ || throw(DimensionMismatch("number of data features must match"))
    k > 1 || throw(ArgumentError("number of clusters must at least be 2"))
    k <= n || throw(ArgumentError("more clusters than data instances are not allowed"))

    C = Matrix{Int}(undef, 0, 0)
    W = Matrix{Float64}(undef, 0, 0)
    Y = zeros(Float64, k, n)
    μ2 = convert(Matrix{Float64}, μ)

    pcs = Vector{PartitionalClustering}(undef, 0)

    distances = pairwise(dist, μ2, X, dims=2)
    pre_objcosts = 0
    objcosts = 0
    for i = 1:maxiter
        # Update cluster assignments, objective function costs, and cluster
        # centers (part 1/2)
        @inbounds for j = 1:n
            cost, y = findmin(view(distances, :, j))
            Y[y, j] = 1
            objcosts += cost
            μ2[:, y] += X[:, j]
        end

        # Update cluster centers (part 2/2)
        @inbounds for j = 1:k
            sz = sum(Y[j, :])
            μ2[:, j] = iszero(sz) ? μ[:, j] : μ2[:, j] / sz
        end

        pairwise!(distances, dist, μ2, X, dims=2)

        pc = PartitionalClustering(copy(r), copy(c), copy(C), copy(W), copy(Y),
                (μ=copy(μ2),))
        push!(pcs, pc)

        # Check for convergence
        isapprox(objcosts, pre_objcosts, atol=ϵ) && break

        fill!(Y, zero(eltype(Y)))
        fill!(μ2, zero(eltype(μ2)))

        pre_objcosts = objcosts
        objcosts = 0
    end

    return pcs
end

"""
    kmeans(X::AbstractMatrix{<:Real}, μ::AbstractMatrix{<:Real}; <keyword arguments>)

Like [`kmeans`](@ref), but automatically generate `r` and `c` according to the
size of `X`.
"""
function kmeans(X::AbstractMatrix{<:Real}, μ::AbstractMatrix{<:Real};
                maxiter::Int=256, dist::SemiMetric=SqEuclidean(),
                ϵ::AbstractFloat=1.0e-6)
    r = collect(1:size(X, 1))
    c = collect(1:size(X, 2))
    return kmeans(X, r, c, μ, maxiter=maxiter, dist=dist, ϵ=ϵ)
end


text:
#===================================================================================================
Mutation operators in differential evolution.
===================================================================================================#

"""
    mutate_rand_1(i::Integer, de::StdDE)

DE/rand/1 mutation, which produces a donor vector to crossover with the `i`-th target vector (parent).
"""
function mutate_rand_1(i::Integer, de::StdDE)
    V = de.pop[:, 1]  # copy
    mutate_rand_1!(V, i, de)
    return V
end

"""
    mutate_rand_1!(V::AV{T}, i::Integer, de::StdDE{<:AM{T}}) where T<:AF

DE/rand/1 mutation in place, which produces a donor vector `V` to crossover with the `i`-th target vector (parent).
"""
function mutate_rand_1!(V::AV{T}, i::Integer, de::StdDE{<:AM{T}}) where T<:AF
    mutate_rand_1!(V, i, de.pop, de.options.F, de._rand_indices)
end

"""
    mutate_rand_1!(V::AV{T}, i::Integer, pop::AM{T}, F::AF, indices::AV{<:Integer}) where T<:AF

A lower-level implementation. `V` is the donor vector to be filled, `i` is the index of the target
vector in the population `pop` whose size is ``n``, `F` is the scale factor, and `indices` is an
array of length ``n`` that contains intergers from 1 to ``n`` (in an arbitrary order).
"""
function mutate_rand_1!(V::AV{T}, i::Integer, pop::AM{T}, F::AF, indices::AV{<:Integer}) where T<:AF
    shuffle!(indices)
    # choose three vectors different from `i` for differential operation
    chosen = MVector{3, eltype(indices)}(0, 0, 0)
    count = 0
    for idx in indices
        if idx != i
            count += 1
            chosen[count] = idx
            if count == 3
                break
            end
        end
    end
    r1, r2, r3 = chosen
    @views V .= pop[:, r1] .+ F.*(pop[:, r2] .- pop[:, r3])

end


text:
function instance_check_unit_tests()

    check_if_supports_tests()
    check_positive_range_tests()
    check_time_windows_tests()
    check_repeated_ids_tests()

end

function check_repeated_ids_tests()
    data = RVRP.RvrpInstance()
    data.locations = [RVRP.Location(;id = "repeated_id"),
                      RVRP.Location(;id = "repeated_id")]
    @test_throws ErrorException RVRP.check_repeated_ids(data)
end

function check_if_supports_tests()

    supported_vec = [BitSet([1, 4, 10, 20, 30, 100]), BitSet([1, 2, 3, 4, 10, 20, 30, 100])]
    features = BitSet([1, 2, 10, 20])
    @test RVRP.check_if_supports(supported_vec, features) == true

    features = BitSet([1, 2, 23])
    @test RVRP.check_if_supports(supported_vec, features) == false

end

function check_positive_range_tests()
    range = RVRP.Range()
    @test RVRP.check_positive_range(range, "") == nothing

    range = RVRP.Range(-1.0, 10.0)
    @test_throws ErrorException RVRP.check_positive_range(range, "")

    range = RVRP.Range(1.0, -10.0)
    @test_throws ErrorException RVRP.check_positive_range(range, "")

    range = RVRP.Range(100.0, 10.0)
    @test_throws ErrorException RVRP.check_positive_range(range, "")
end

function check_time_windows_tests()

    tws = [RVRP.Range(10.0, 20.0), RVRP.Range(15.0, 20)]
    @test_throws ErrorException RVRP.check_time_windows(tws, "")

    tws = [RVRP.Range(10.0, 20.0), RVRP.Range(0.0, 5.0)]
    @test_throws ErrorException RVRP.check_time_windows(tws, "")

end


text:
using FactCheck, MarketData
FactCheck.setstyle(:compact)
FactCheck.onlystats(true)

facts("first methods takes the first observations") do

    context("defaults to n=1") do
        @fact first(cl).values[1] --> 111.94
    end

    context("takes other n values") do
        @fact first(cl, 2).values[2] --> 102.5
    end

end


text:
using Homotopies
using Base.Test
import FixedPolynomials
const FP = FixedPolynomials
import DynamicPolynomials: @polyvar

include("straightline_test.jl")
include("geodesic_test.jl")
include("misc_test.jl")
include("condition_test.jl")


text:
using LowRank

using JLD, Test

# rank 1
file = load(joinpath(@__DIR__, "..","datasets", "uniform_2d_10_pts.jld"))
data_set, (err, indices) = file["data"], file["size4_rank1"]
k = 1
n = 4

# @test find_low_rank_subset_checkall(data_set,n,k)[1] <= err + 0.0001
@test find_low_rank_subset_sample_rep(data_set,n,k)[1] <= 2*err + .0001
@test FastLowRank(data_set, n, k, verbose = true)[1] <= err + 0.0001

FastLowRank(data_set, n, k, verbose = true, num_trajectories = 100, convergence_threshold = 0.00000001)

# PLotting
# using Plots
# p = scatter([Tuple(x) for x in data_set], label = "all points",  series_annotations = text.(1:10, :bottom))
# scatter!(p, [Tuple(data_set[i]) for i in indices], label = "optimal points")


text:
#=----------------------------------------------------------------------------

Perceptually Uniform Colour Maps

Copyright (c) 2015 Peter Kovesi
pk@peterkovesi.com

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, subject to the following conditions:

The above copyright notice and this permission notice shall be included in
all copies or substantial portions of the Software.

The Software is provided "as is", without warranty of any kind.

PK October 2015

----------------------------------------------------------------------------=#
"""
**PerceptualColourMaps**

Perceptually Uniform Colour maps

Peter Kovesi

[peterkovesi.com](http://peterkovesi.com)

Reference:
[Good Colour Maps: How to Design Them.  arXiv:1509.03700 [cs.GR] 2015.](http://arxiv.org/abs/1509.03700)

"""
module PerceptualColourMaps

include("cmap.jl")
include("applycolourmaps.jl")
include("viewlabspace.jl")
include("relief.jl")
include("utilities.jl")

end  # module


text:
# StreetRaptor combines RAPTOR with a street search from a geographic origin to a geographic destination

using .OSRM
using ProgressBars
using Dates

struct EgressTime
    dest_idx::Int64
    time_seconds::Int32
end

struct EgressTimes
    destinations::Vector{Coordinate}
    stop_dest_time::Vector{Vector{EgressTime}}
end

struct StreetRaptorRequest
    origin::Coordinate
    departure_time::Int32
    date::Date
    max_access_distance_meters::Float64
    walk_speed_meters_per_second::Float64
    max_rides::Int64
end

struct StreetRaptorResult
    times_at_destinations::Vector{Int32}
    egress_stop_for_destination::Vector{Int64}
    raptor_result::RaptorResult
    request::StreetRaptorRequest
end

# find the egress times from all stops to all destinations, within distance_limit
function find_egress_times(net::TransitNetwork, osrm::OSRMInstance, destinations::Vector{Coordinate}, max_distance_meters::Real)::EgressTimes
    stop_dest_time = Vector{Vector{EgressTime}}()
    sizehint!(stop_dest_time, length(net.stops))
    lat_diff = meters_to_degrees_lat(max_distance_meters)
    for (i, stop) in ProgressBar(enumerate(net.stops))
        origin = Coordinate(stop.stop_lat, stop.stop_lon)
        # prefilter to just destinations that are nearby
        candidate_destinations = bbox_filter(origin, destinations, max_distance_meters)
        egress_times = Vector{EgressTime}()

        if length(candidate_destinations) > 0
            destination_coords = destinations[candidate_destinations]
            dists = distance_matrix(osrm, [origin], destination_coords)

            for candidate_dest_index in 1:length(candidate_destinations)
                time = dists.durations[1, candidate_dest_index]
                dist = dists.distances[1, candidate_dest_index]
                dest_idx = candidate_destinations[candidate_dest_index]
                if dist <= max_distance_meters
                    push!(egress_times, EgressTime(dest_idx, convert(Int32, round(time))))
                end
            end
        end

        push!(stop_dest_time, egress_times)
    end

    @assert length(stop_dest_time) == length(net.stops)

    return EgressTimes(destinations, stop_dest_time)
end

function street_raptor(net::TransitNetwork, access_router::OSRMInstance, req::StreetRaptorRequest, destinations::EgressTimes)::StreetRaptorResult
    @info "performing access search"

    # find stops near origin
    stop_coords = map(s -> Coordinate(s.stop_lat, s.stop_lon), net.stops)
    stops_near_origin = bbox_filter(req.origin, stop_coords, req.max_access_distance_meters)

    access = distance_matrix(access_router, [req.origin], stop_coords[stops_near_origin])

    accessible_stops = Vector{StopAndTime}()
    for stop_near_origin_idx in 1:length(stops_near_origin)
        stop_idx = stops_near_origin[stop_near_origin_idx]
        # this is the time the stop is reached
        time = req.departure_time + access.durations[1, stop_near_origin_idx]
        dist = access.distances[1, stop_near_origin_idx]

        if dist <= req.max_access_distance_meters
            push!(accessible_stops, StopAndTime(stop_idx, round(time)))
        end
    end

    @info "$(length(accessible_stops)) stops found near origin"
    @info "begin transit routing"

    rreq = RaptorRequest(
        accessible_stops,
        req.max_rides,
        req.date,
        req.walk_speed_meters_per_second
    )

    raptor_res = raptor(net, rreq)

    @info "transit routing complete. adding egress times."
    times_at_destinations::Vector{Int32} = fill(MAX_TIME, length(destinations.destinations))
    egress_stops::Vector{Int64} = fill(INT_MISSING, length(destinations.destinations))

    for stopidx in 1:length(net.stops)
        time_at_stop = raptor_res.times_at_stops_each_round[size(raptor_res.times_at_stops_each_round, 1), stopidx]
        if time_at_stop < MAX_TIME
            # this stop was reached
            for egress in destinations.stop_dest_time[stopidx]
                time_at_dest = time_at_stop + egress.time_seconds
                if time_at_dest < times_at_destinations[egress.dest_idx]
                    times_at_destinations[egress.dest_idx] = time_at_dest
                    egress_stops[egress.dest_idx] = stopidx
                end
            end
        end
    end

    return StreetRaptorResult(times_at_destinations, egress_stops, raptor_res, req)
end


text:
import Base: collect, getindex, length, maximum, minimum, size, NamedTuple

AbstractSample{DT} = Union{Nothing, AbstractVector{DT}}

_sort(s::Nothing) = s
_sort(s::AbstractVector) = sort(unique(collect(s)))

"""

"""
struct Parameter{DT <: Number, ST <: AbstractSample{DT}}
    name::Symbol
    minimum::DT
    maximum::DT
    samples::ST

    function Parameter(name::Symbol, minimum::DT, maximum::DT, samples::ST) where {DT, ST <: AbstractSample{DT}}
        @assert minimum ≤ maximum
        if typeof(samples) <: AbstractVector
            @assert all([minimum ≤ s for s in samples])
            @assert all([maximum ≥ s for s in samples])
        end
        _samples = _sort(samples)
        new{DT, typeof(_samples)}(name, minimum, maximum, _samples)
    end
end

Parameter(name::AbstractString, minimum, maximum, samples) = Parameter(Symbol(name), minimum, maximum, samples)
Parameter(name, minimum::DT, maximum::DT) where {DT} = Parameter(name, minimum, maximum, nothing)
Parameter(name, minimum::DT, maximum::DT, n::Int) where {DT} = Parameter(name, minimum, maximum, LinRange(minimum, maximum, n))
Parameter(name, samples::AbstractVector) = Parameter(name, minimum(samples), maximum(samples), samples)

Base.hash(p::Parameter, h::UInt) = hash(p.name, hash(p.minimum, hash(p.maximum, hash(p.samples, h))))

Base.:(==)(p1::Parameter, p2::Parameter) = (
                                p1.name    == p2.name
                             && p1.minimum == p2.minimum
                             && p1.maximum == p2.maximum
                             && p1.samples == p2.samples)

function show(io::IO, p::Parameter)
    println(io, "Parameter $(p.name) with ")
    println(io, "   minimum = ", p.minimum)
    println(io, "   maximum = ", p.maximum)
    println(io, "   samples = ")
    show(io, p.samples)
end

const ParameterWithSamples = Parameter{DT,ST} where {DT, ST <: AbstractVector}
const ParameterWithoutSamples = Parameter{DT,ST} where {DT, ST <: Nothing}

Base.collect(p::Parameter) = p.samples
Base.maximum(p::Parameter) = p.maximum
Base.minimum(p::Parameter) = p.minimum

Base.length(p::ParameterWithSamples) = length(p.samples)
Base.length(p::ParameterWithoutSamples) = 0

Base.size(p::ParameterWithSamples) = size(p.samples)
Base.size(p::ParameterWithoutSamples) = (0,)

hassamples(p::ParameterWithSamples) = length(p.samples) > 0
hassamples(p::ParameterWithoutSamples) = false

@inline Base.@propagate_inbounds Base.getindex(p::ParameterWithSamples, i) = p.samples[i]
@inline Base.@propagate_inbounds Base.getindex(p::ParameterWithoutSamples, i) = error("Parameter $(p.name) indexed with $(i) but has no samples.")
@inline Base.@propagate_inbounds Base.getindex(p::ParameterWithSamples, ::Colon) = p.samples
@inline Base.@propagate_inbounds Base.getindex(p::ParameterWithoutSamples, ::Colon) = error("Parameter $(p.name) has no samples.")

function Base.NamedTuple(parameters::Vararg{Parameter{DT}}) where {DT}
    names = Tuple(p.name for p in parameters)
    NamedTuple{names}(parameters)
end

function Parameter(h5::H5DataStore, path::AbstractString)
    g = h5[path]
    name = _name(g)

    minimum = read(g["minimum"])
    maximum = read(g["maximum"])
    samples = read(g["samples"])

    Parameter(name, minimum, maximum, samples)
end

function h5save(param::Parameter, h5::H5DataStore, path::AbstractString=string(param.name))
    g = _create_group(h5, path)
    g["minimum"] = param.minimum
    g["maximum"] = param.maximum
    g["samples"] = param.samples
end

function h5load(::Type{Parameter}, h5::H5DataStore, path::AbstractString)
    Parameter(h5, path)
end


text:
include("../src/expmatrices.jl")
using LinearAlgebra
using Random
using Test
function testexp()
    A=zeros(Int64, 4, 4)
    A[1,3] = 1
    A[3,1] = -1
    @time @testset "test 1 exp big matrices" begin
        @test isapprox( one(A), _expmat1(0A), atol=1e-16)
        @test isapprox(exp(0.123*A), _expmat1(0.123*A), atol=1e-15)
    end
    @time @testset "test 2 exp big matrices" begin
        Abig = convert(Array{BigFloat,2},A)
        Random.seed!(81725)
        sens=1
        for k=1:4
            setprecision(k*256) do
                for i = 1:100:10000
                    v=rand(BigFloat)*i*sens
                    res = _expmat0(v*Abig)
                    resRef = one(Abig)
                    resRef[1,1] = resRef[3,3] = cos(v)
                    resRef[1,3] = sin(v)
                    resRef[3,1]= -resRef[1,3]
        #           println("i=$i norm=$(norm(resRef-res))")
                    @test isapprox( resRef, res, atol=(1e-77)^k)
                    sens *= -1
                end
            end
        end
    end
    @time @testset "test 3 exp big matrices" begin
        Random.seed!(8781115)

        b = [ 0 0 0 1; 0 0 -1 0; 0 1 0 0; -1 0 0 0]
        for i=1:20
            C=rand(4,4)*i
            B = b*i
            @test isapprox(exp(C), _expmat1(C), rtol=1e-13)
            @test isapprox(exp(0.7C), _expmat1(0.7C), rtol=1e-12)
            @test isapprox(exp(0.89C), _expmat1(0.89C), rtol=1e-12)
            @test isapprox(exp(B), _expmat1(B), rtol=1e-14)
            @test isapprox(exp(10B), _expmat1(10B), rtol=1e-13)
            @test isapprox(exp(100B), _expmat1(100B), rtol=1e-12)
            @test isapprox(exp(1000B), _expmat1(1000B), rtol=1e-11)
            @test isapprox(exp(10000B), _expmat1(10000B), rtol=1e-10)
        end
    end
    @time @testset "test 4 exp big matrices big precision" begin
        Random.seed!(90908)
        prec_old = prec = precision(BigFloat)
        for i=1:6
            prec *= 2
            setprecision(prec)
            @test isapprox(_expmat1(2big(pi)*A)-I,zeros(BigFloat,4,4),atol=eps(BigFloat)*100)
        end
        setprecision(prec_old)
    end
end
testexp()


text:
using OpenQuantumBase, Test

A = (s) -> (1 - s)
B = (s) -> s
u = [1.0 + 0.0im, 1] / sqrt(2)
ρ = u * u'

H = DenseHamiltonian([A, B], [σx, σz])

@test_throws ArgumentError DenseHamiltonian([A, B], [σx, σz], unit=:hh)
@test size(H) == (2, 2)
@test H(0) == 2π * σx
@test evaluate(H, 0) == σx
@test H(0.5) == π * (σx + σz)
@test evaluate(H, 0.5) == (σx + σz) / 2
@test get_cache(H) ≈ π * (σx + σz)

# update_cache method
C = similar(σz)
update_cache!(C, H, 10, 0.5)
@test C == -1im * π * (σx + σz)

# update_vectorized_cache method
C = get_cache(H)
C = C⊗C
update_vectorized_cache!(C, H, 10, 0.5)
temp = -1im * π * (σx + σz)
@test C == σi ⊗ temp - transpose(temp) ⊗ σi

# in-place update for matrices
du = [1.0 + 0.0im 0; 0 0]
H(du, ρ, 2, 0.5)
@test du ≈ -1.0im * π * ((σx + σz) * ρ - ρ * (σx + σz))

# eigen-decomposition
w, v = eigen_decomp(H, 0.5)
@test w ≈ [-1, 1] / sqrt(2)
w, v = eigen_decomp(H, 0.0)
@test w ≈ [-1, 1]
@test v ≈ [-1 1; 1 1] / sqrt(2)