using Plots
gr(); ENV["PLOTS_TEST"] = "true"
#clibrary(:colorcet)
clibrary(:misc)

function pngplot(P...; kwargs...)
    sleep(0.1)
    pngfile = tempname() * ".png"
    savefig(plot(P...; kwargs...), pngfile)
    showimg("image/png", pngfile)
end
pngplot(; kwargs...) = pngplot(plot!(; kwargs...))

showimg(mime, fn) = open(fn) do f
    base64 = base64encode(f)
    display("text/html", """<img src="data:$mime;base64,$base64">""")
end

using SymPy
#sympy[:init_printing](order="lex") # default
#sympy[:init_printing](order="rev-lex")

using SpecialFunctions
using QuadGK

D(q,p) = q*log(q/p) + (1-q)*log((1-q)/(1-p))
q = 0:0.01:1
PP = []
for p in [0.2, 0.4, 0.6, 0.8]
    P = plot(legend=false)
    plot!(title="y = D(q, $(round(p,3)))", titlefontsize=10)
    plot!(xlabel="q_1", ylabel="y")
    plot!(q, D.(q, p))
    push!(PP, P)
end
plot(PP..., size=(700, 500))

D(q,p) = q*log(q/p) + (1-q)*log((1-q)/(1-p))
p = 0:0.01:1
PP = []
for q in [0.2, 0.4, 0.6, 0.8]
    P = plot(legend=false)
    plot!(title="y = D($(round(q,3)), p)", titlefontsize=10)
    plot!(xlabel="p_1", ylabel="y")
    plot!(p, D.(q, p))
    push!(PP, P)
end
plot(PP..., size=(700, 500))

pdflogistic(μ,s,x) = 1/(4s)*sech((x-μ)/(2s))^2
pdfnormal(μ,σ,x) = exp(-(x-μ)^2/(2*σ^2))/√(2π*σ^2)

μ = 0.0
s = 1.0
σ = π*s/√3
x = -10:0.02:10

P = plot(legend=:topleft, ylims=(0, 0.26))
plot!(x, pdflogistic.(μ,s,x), label="logistic", lw=2)
plot!(x, pdfnormal.(μ,σ,x),   label="normal",    lw=2, ls=:dash)

binomcoeff(n,k) = exp(lgamma(n+1)-lgamma(k+1)-lgamma(n-k+1))
pdfbinom(n,p,k) = 0 ≤ k ≤ n ? binomcoeff(n,k)*p^k*(1-p)^(n-k) : 0.0
pdfnormal(n,p,x) = exp(-x^2/(2p*(1-p)))/√(2π*n*p*(1-p))

PP = []
for (n,p) in [(10,0.4), (25,0.4), (50,0.4), (100,0.4)]
    k = 0:n
    kk = 0:n/400:n
    x = (kk .- n*p)./√n
    P = plot(legend=:topright, xlabel="k")
    plot!(k, pdfbinom.(n,p,k),  label="Bin($n, $p)",  lw=2)
    plot!(kk, pdfnormal.(n,p,x), label="Normal dist", lw=2, ls=:dash)
    push!(PP, P)
end
plot(PP..., size=(800, 600))
#pngplot(PP..., size=(800, 600))