using Base.MathConstants
using Base64
using Printf
using Statistics
const e = ℯ
endof(a) = lastindex(a)
linspace(start, stop, length) = range(start, stop, length=length)

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

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
using SpecialFunctions
using QuadGK

# Override the Base.show definition of SymPy.jl:
# https://github.com/JuliaPy/SymPy.jl/blob/29c5bfd1d10ac53014fa7fef468bc8deccadc2fc/src/types.jl#L87-L105

@eval SymPy function Base.show(io::IO, ::MIME"text/latex", x::SymbolicObject)
    print(io, as_markdown("\\displaystyle " * sympy.latex(x, mode="plain", fold_short_frac=false)))
end
@eval SymPy function Base.show(io::IO, ::MIME"text/latex", x::AbstractArray{Sym})
    function toeqnarray(x::Vector{Sym})
        a = join(["\\displaystyle " * sympy.latex(x[i]) for i in 1:length(x)], "\\\\")
        """\\left[ \\begin{array}{r}$a\\end{array} \\right]"""
    end
    function toeqnarray(x::AbstractArray{Sym,2})
        sz = size(x)
        a = join([join("\\displaystyle " .* map(sympy.latex, x[i,:]), "&") for i in 1:sz[1]], "\\\\")
        "\\left[ \\begin{array}{" * repeat("r",sz[2]) * "}" * a * "\\end{array}\\right]"
    end
    print(io, as_markdown(toeqnarray(x)))
end

function bisection(f, a, b; atol=eps(), maxiter=10^3)
    f(a)*f(b) > 0 && return NaN
    aₙ, bₙ = a, b
    k = 0
    while abs(aₙ - bₙ) > atol
        k > maxiter && return NaN
        c = (aₙ+bₙ)/2
        if f(bₙ)*f(c) ≥ 0
            aₙ, bₙ = aₙ, c
        else
            aₙ, bₙ = c, bₙ
        end
        k += 1
    end
    return (aₙ+bₙ)/2
end

f(x) = x^2 - 2
@show ξ = bisection(f, 1, 2)
@show √2
@show round(f(ξ), digits=15)

x = 1:0.005:2
plot(size=(400,250), legend=:topleft)
plot!(x, f.(x), label="y = x^2 - 2")
hline!([0], label="y = 0")
vline!([ξ], label="x = $(round(ξ, digits=5))", ls=:dash)

f(x) = cos(x) - x
@show ξ = bisection(f, 0, 1)
@show round(f(ξ), digits=15)

x = 0:0.005:1
plot(size=(400,250))
plot!(x, f.(x), label="y = cos x - x")
hline!([0], label="y = 0")
vline!([ξ], label="x = $(round(ξ, digits=5))", ls=:dash)

f(x) = 0.6 + 0.4x - 0.8x^2 + 0.3sin(10x-0.5)
ξ = bisection(x->f(x)-x, 0, 1)
x = 0:0.01:1
plot(size=(400,400), aspect_ration=1, legend=false)
plot!(xlim=(0,1), ylim=(0,1))
plot!(x, f.(x), lw=1.5)
plot!([0,1], [0,1], color=:black)
hline!([ξ], ls=:dash, color=:red)
vline!([ξ], ls=:dash, color=:red)

φ(x) = sin(x) - 0.5x + 0.1x^2
x = -2.5:0.02:3.5
ε = 0.3
plot(size=(500, 350))
plot!(x, φ.(x), label="y = phi(x)", color=:black, lw=2)
plot!(x, φ.(x).+ε, label="phi(x) - eps < y < phi(x)+eps", color=:red, ls=:dash, fill=(:pink, 0.2, φ.(x)))
plot!(x, φ.(x).-ε, label="", color=:red, ls=:dash, fill=(:pink, 0.2, φ.(x)))

f(x) = sin(x) - 0.5x + 0.1x^2
g(x) = cos(x) - 0.5 + 0.1x^3
plot(size=(450, 280), legend=:top)
x = -2.5:0.02:3.5
plot!(x, f.(x), color=:blue, label="f(x)")
plot!(x, g.(x), color=:red,  label="g(x)")
x = -2.0:0.02:3.0
plot!(x, g.(x), color=:red, fill=(0.1, f.(x)), label="")
vline!([minimum(x)], color=:black, ls=:dash,    label="x=a")
vline!([maximum(x)], color=:black, ls=:dashdot, label="x=b")

x = symbols("x")
diffexpneginvx(n) = diff(exp(-1/x), x, n)
[diffexpneginvx(n) for n in 0:3]

f(x) = x > 0 ? exp(-1/x) : zero(x)
g(x) = x > 0 ? exp(-1/x)/x^2 : zero(x)

p1 = plot(ylims=(-0.1, 1.0), legend=:topleft)
plot!(f, -10, 20, label="f(x)")

p2 = plot(ylims=(-0.1, 0.6), legend=:topright)
plot!(g, -10, 20, label="df(x)/dx")

plot(p1, p2, size=(700,200))

f(x) = 0 < x < 1 ? exp(-1/x+1/(x-1)) : zero(x)
@show Z, err = quadgk(f, 0, 1) # = ∫_0^1 f(x) dx, error
g(x) = f(x)/Z
h(x) = quadgk(g, -1, x)[1] # = ∫_{-1}^x g(t) dt

p1 = plot(legend=:topleft, ylims=(-0.1,2.7))
plot!(g, -1, 2, label="g(x)")

p2 = plot(legend=:topleft, ylims=(-0.1, 1.1))
plot!(h, -1, 2, label="h(x)")

plot(p1, p2, size=(700,200))

φ(x) = h(x)*h(4-x)

plot(size=(400, 200), ylims=(-0.1, 1.1))
plot!(φ, -1, 5, label="phi(x)")

f(x) = x == 0 ? 0 : sin(1/x)

p1 = plot(size=(400, 250), xlabel="x", ylabel="y")
plot!(f, -1.0, 1.0, ylims=(-1.2, 2.0), label="y = sin(1/x)")

p2 = plot(size=(400, 250), xlabel="x", ylabel="y")
plot!(f, -0.1, 0.1, ylims=(-1.2, 2.0), label="y = sin(1/x)")

plot(p1, p2, size=(700, 250))

f(x) = x == 0 ? 0 : x*sin(1/x)

p1 = plot(size=(400, 250), xlabel="x", ylabel="y")
plot!(f, -3, 3, ylims=(-0.4, 1.5), label="y = x sin(1/x)")

p2 = plot(size=(400, 250), xlabel="x", ylabel="y")
plot!(f, -0.3, 0.3, ylims=(-0.25, 0.3), label="y = x sin(1/x)")

plot(p1, p2, size=(700, 250))

struct WeierstrassFunction{R<:Real, S<:Integer, T<:Integer} <: Function
    a::R;  b::S;  N::T
end
(F::WeierstrassFunction)(x) = sum(n->F.a^n*cos(F.b^n*π*x), 0:F.N)

F = WeierstrassFunction(0.6, 13, 100)
@show F.a*F.b, 1+3π/2

p1 = plot(legend=false, titlefontsize=11)
plot!(title="Weierstrass function: a = $(F.a), b = $(F.b)")
plot!(size=(600, 420), xlabel="x", ylabel="y")
plot!(F, -2, 2)

p2 = plot(legend=false, titlefontsize=11)
plot!(title="Weierstrass function: a = $(F.a), b = $(F.b)")
plot!(size=(600, 420), xlabel="x", ylabel="y")
plot!(F, -2/F.b, 2/F.b)

plot(p1, p2, size=(700, 250))

F = WeierstrassFunction(0.85, 7, 100)
@show F.a*F.b, 1+3π/2

p1 = plot(legend=false, titlefontsize=11)
plot!(title="Weierstrass function: a = $(F.a), b = $(F.b)")
plot!(size=(600, 420), xlabel="x", ylabel="y")
plot!(F, -2, 2)

p2 = plot(legend=false, titlefontsize=11)
plot!(title="Weierstrass function: a = $(F.a), b = $(F.b)")
plot!(size=(600, 420), xlabel="x", ylabel="y")
plot!(F, -2/F.b, 2/F.b)

plot(p1, p2, size=(700, 250))

F = WeierstrassFunction(1/4, 3, 100)
@show F.a*F.b, 1+3π/2

p1 = plot(legend=false, titlefontsize=11)
plot!(title="Weierstrass function: a = $(F.a), b = $(F.b)")
plot!(size=(600, 420), xlabel="x", ylabel="y")
plot!(F, -2, 2)

p2 = plot(legend=false, titlefontsize=11)
plot!(title="Weierstrass function: a = $(F.a), b = $(F.b)")
plot!(size=(600, 420), xlabel="x", ylabel="y")
plot!(F, -2/F.b, 2/F.b)

plot(p1, p2, size=(700, 250))

seesaw(x) = abs(x - round(x))
plot(seesaw, size=(400, 250), ylims=(-0.1, 0.7), label="seesaw(x)")

struct TakagiFunction{R<:Real, T<:Integer} <: Function
    w::R;  N::T
end
(F::TakagiFunction)(x) = sum(n->F.w^n*seesaw(2^n*x), 0:F.N)

F = TakagiFunction(0.5, 100)

p1 = plot(legend=false, titlefontsize=11)
plot!(title="Takagi function: w = $(F.w)")
plot!(size=(600, 420), xlabel="x", ylabel="y")
x = 0:0.001:1
plot!(x, F.(x))

p2 = plot(legend=false, titlefontsize=11)
plot!(title="Takagi function: w = $(F.w)")
plot!(size=(600, 420), xlabel="x", ylabel="y")
x = 0.25:0.00025:0.50
plot!(x, F.(x))

plot(p1, p2, size=(700, 250))

F = TakagiFunction(0.7, 100)

p1 = plot(legend=false, titlefontsize=11)
plot!(title="Takagi function: w = $(F.w)")
plot!(size=(600, 420), xlabel="x", ylabel="y")
x = 0:0.001:1
plot!(x, F.(x))

p2 = plot(legend=false, titlefontsize=11)
plot!(title="Takagi function: w = $(F.w)")
plot!(size=(600, 420), xlabel="x", ylabel="y")
x = 0.250:0.00025:0.50
plot!(x, F.(x))

plot(p1, p2, size=(700, 250))

F = TakagiFunction(1/4, 100)

p1 = plot(legend=false, titlefontsize=11)
plot!(title="Takagi function: w = $(F.w)")
plot!(size=(600, 420), xlabel="x", ylabel="y")
x = 0:0.001:1
plot!(x, F.(x))

p2 = plot(legend=false, titlefontsize=11)
plot!(title="y = x(1-x)")
plot!(size=(600, 420), xlabel="x", ylabel="y")
x = 0:0.001:1
plot!(x, @.(x*(1-x)))

plot(p1, p2, size=(700, 250))

f(x) = 1/x
x = 0.1:0.001:1

ε = 0.7

plot(xlims=(0,1), ylims=(0,10), size=(500, 350))
plot!(x, f.(x), label="y = 1/x", lw=2)

a = 0.15
δ₀ = a - 1/(1/a+ε)
δ₁ = 1/(1/a-ε) - a
c = :red
plot!([a-δ₀,a-δ₀], [0,f(a-δ₀)], color=c, label="")
plot!([a,a], [0,f(a)], color=c, lw=2, label="x = $a")
plot!([a+δ₁,a+δ₁], [0,f(a+δ₁)], color=c, label="")
plot!([0,a-δ₀], [f(a-δ₀),f(a-δ₀)], ls=:dash, color=c, label="y = 1/$a+$ε")
plot!([0,a], [f(a),f(a)], color=c, ls=:dash, lw=2, label="y = 1/$a")
plot!([0,a+δ₁], [f(a+δ₁),f(a+δ₁)], ls=:dashdot, color=c, label="y = 1/$a-$ε")

a = 0.5
δ₀ = a - 1/(1/a+ε)
δ₁ = 1/(1/a-ε) - a
c = :orange
plot!([a-δ₀,a-δ₀], [0,f(a-δ₀)], color=c, label="")
plot!([a,a], [0,f(a)], color=c, lw=2, label="x = $a")
plot!([a+δ₁,a+δ₁], [0,f(a+δ₁)], color=c, label="")
plot!([0,a-δ₀], [f(a-δ₀),f(a-δ₀)], ls=:dash, color=c, label="y = 1/$a+$ε")
plot!([0,a], [f(a),f(a)], color=c, ls=:dash, lw=2, label="y = 1/$a")
plot!([0,a+δ₁], [f(a+δ₁),f(a+δ₁)], ls=:dashdot, color=c, label="y = 1/$a-$ε")

f(n,x) = x^n

x = 0:0.005:1
p1 = plot(xlims=(0,1), legend=:topleft)
cycls = [:solid, :dash, :dashdot, :dashdotdot]
for n in 1:12
    plot!(x, f.(n,x), label="n = $n", ls=cycls[mod1(n, lastindex(cycls))])
end
plot(p1, size=(360,370), aspect_ratio=1, title="y = x^n", lw=1.5)

f(n,x) = x^n

x = 0:0.005:1

n = 3
p1 = plot(xlims=(-0.01,1.01), legend=:topleft, aspect_ratio=1)
plot!(x, f.(n,x), label="y = x^$n")
plot!(x[1:end-1], 1/2*fill(1, size(x[1:end-1])), label="", color=:red, ls=:dash, fill=(0.1, 0))
plot!([1,1], [1/2,1], label="", color=:red, alpha=0.1, lw=3)
plot!([1/2^(1/n), 1/2^(1/n)], [0, 1/2], label="x = 1/2^{1/$n}", color=:black, ls=:dash)

n = 9
p2 = plot(xlims=(-0.01,1.01), legend=:topleft, aspect_ratio=1)
plot!(x, f.(n,x), label="y = x^$n")
plot!(x[1:end-1], 1/2*fill(1, size(x[1:end-1])), label="", color=:red, ls=:dash, fill=(0.1, 0))
plot!([1,1], [1/2,1], label="", color=:red, alpha=0.1, lw=3)
plot!([1/2^(1/n), 1/2^(1/n)], [0, 1/2], label="x = 1/2^{1/$n}", color=:black, ls=:dash)

plot(p1, p2, size=(600, 250))

g(N, s) = sum(n->1/n^s, 1:N)
s = 1.01:0.01:5.00

p1 = plot(title="sum of 1/n^s for n=1,...,N")
plot!(ylims=(0.8, 5), xlabel="s")
plot!(s, zeta.(s), label="zeta(s)", lw=2)
for (ls, N) in [(:dot,3), (:dash, 10), (:dashdot, 30), (:dashdotdot, 100), (:dashdotdot, 300), (:dashdotdot, 1000)]
    plot!(s, g.(N,s), label="N = $N", ls=ls, lw=1.5)
end
plot(p1, size=(500, 350))

peano_x₀(t) = t < 1/2 ? 1/2 : t
peano_y₀(t) = t < 1/2 ? t : 1/2

n = 0
t = 0:1/(2*4^n):1
x = peano_x₀.(t)
y = peano_y₀.(t)
plot(size=(200,200))
plot!(title="Peano($n)", titlefontsize=10)
plot!(xlims=(0,1), ylims=(0,1))
plot!(aspect_ratio=1, legend=:false)
plot!(x, y)

function peano_x(n, t)
    if n == 0
        peano_x₀(t)
    else
        if t < 1/(2*4^n)
            1/(2*2^n)
        elseif t < 1/4
            peano_y(n-1, 4t)/2
        elseif t < 2/4
            peano_x(n-1, 4t-1)/2
        elseif t < 3/4
            (1-peano_x(n-1, 3-4t))/2 + 1/2
        else
            (1-peano_y(n-1, 4-4t))/2 + 1/2
        end
    end
end

function peano_y(n, t)
    if n == 0
        peano_y₀(t)
    else
        if t < 1/(2*4^n)
            2^n*t
        elseif t < 1/4
            peano_x(n-1, 4t)/2
        elseif t < 2/4
            peano_y(n-1, 4t-1)/2 + 1/2
        elseif t < 3/4
            peano_y(n-1, 3-4t)/2 + 1/2
        else
            peano_x(n-1, 4-4t)/2
        end
    end
end

P = Plots.Plot[]
for n in 0:8
    t = 0:1/(2*4^n):1
    x = peano_x.(n,t)
    y = peano_y.(n,t)
    p = plot()
    plot!(title="Peano($n)", titlefontsize=10)
    plot!(xlims=(0,1), ylims=(0,1))
    plot!(aspect_ratio=1, legend=:false)
    plot!(x, y)
    push!(P, p)
end
plot(P..., size=(800,800), layout=@layout(grid(3,3)))
pngplot()

P = Plots.Plot[]
for n in 0:8
    t = 0:1/(2*4^n):1
    x = peano_x.(n,t)
    y = peano_y.(n,t)
    p = plot()
    plot!(title="Peano($n)", titlefontsize=10)
    plot!(xlims=(0,1), ylims=(0,1))
    plot!(aspect_ratio=1, legend=:false)
    plot!(t, x)
    plot!(t, y)
    push!(P, p)
end
plot(P..., size=(800,800), layout=@layout(grid(3,3)))
pngplot()

n = 8
px = plot()
py = plot()
t = 0:1/(2*4^n):1
x = peano_x.(n,t)
y = peano_y.(n,t)
plot!(px, title="Peano_x($n)", titlefontsize=10)
plot!(py, title="Peano_y($n)", titlefontsize=10)
plot!(px, xlims=(0,1), ylims=(0,1))
plot!(py, xlims=(0,1), ylims=(0,1))
plot!(px, legend=:topleft)
plot!(py, legend=:topleft)
plot!(px, t, x, label="$n", lw=0.3)
plot!(py, t, y, label="$n", lw=0.3)
plot(px, py, size=(600,500), layout=@layout(grid(2,1)), legend=false)
pngplot()

P = Plots.Plot[]
for n in 1:9
    t = 0:1/(2*4^n):1
    dx = peano_x.(n,t) - peano_x.(n-1,t)
    dy = peano_y.(n,t) - peano_y.(n-1,t)
    p = plot()
    plot!(title="Peano($n) - Peano($(n-1))", titlefontsize=10)
    plot!(xlims=(0,1), ylims=(-1/2,1/2))
    plot!(legend=:false)
    plot!(t, dx, lw=0.7)
    plot!(t, dy, lw=0.7)
    push!(P, p)
end
plot(P..., size=(800,600))
pngplot()