print("using Plots ->")
@time using Plots
gr(legend=false); ENV["PLOTS_TEST"] = "true"

print("\nplot(sin) ->")
@time plot(sin, size=(300, 160)) |> display # コンパイル

print("\nusing PyPlot ->")
@time using PyPlot: PyPlot, plt

print("\nusing DifferentialEquations ->")
@time using DifferentialEquations

using Base64
using Combinatorics
using Distributed
using Distributions
using Libdl
using LinearAlgebra
using Optim
using Primes
using ProgressMeter
using Random: seed!
using RCall
using SpecialFunctions
using Statistics
using SymPy: SymPy, sympy, Sym, @vars, @syms, simplify, oo, PI

ldisp(x...) = display("text/html", raw"$$" * prod(x) * raw"$$")

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

### 函数にせずに実行

@time begin
    L = 10^7
    c = 0
    for i in 1:L
        global c
        c += ifelse(rand()^2+rand()^2 ≤ 1, 1, 0)
    end
    4c/L
end

### 函数にして実行

function simpi(L=10^7)
    c = 0
    for i in 1:L
        c += ifelse(rand()^2+rand()^2 ≤ 1, 1, 0)
    end
    4c/L
end
simpi(10^5) # simpi(::Int64)がコンパイルされる

@time simpi()

versioninfo()

run(`gcc --version`); # ほとんど使っていないのでバージョンが古い

### gcc

C_code= raw"""
#include <time.h>
#include <stdlib.h>

double simpi(unsigned long n){
    double x,y;
    srand((unsigned)time(NULL));
    double R = (double) RAND_MAX;
    unsigned long count = 0;
    for(unsigned long j=0; j < n; j++){
        x = ((double) rand())/R;
        y = ((double) rand())/R;
        if(x*x + y*y <= 1){
            count++;
        }
    }
    return ((double) 4.0)*((double) count)/((double) n);
}
"""

filename = tempname()
filenamedl = filename * "." * Libdl.dlext

open(`gcc -Wall -O3 -march=native -xc -shared -o $filenamedl -`, "w") do f
     print(f, C_code)
end

## run(`ls -l $filenamedl`);

const simpi_gcc_rand_lib = filename
simpi_gcc_rand(n::Int64) = ccall((:simpi, simpi_gcc_rand_lib), Float64, (Int64,), n)
simpi_gcc_rand(10^5)

@time simpi_gcc_rand(10^8)

### Julia

@time simpi(10^8)

### gcc with dSFMT

## Download http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/SFMT/dSFMT-src-2.2.3.tar.gz
##
## Extract
##   dSFMT-src-2.2.3/dSFMT.h
##   dSFMT-src-2.2.3/dSFMT.c

C_code= raw"""
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#define HAVE_SSE2
#define DSFMT_MEXP 521
#include "dSFMT-src-2.2.3/dSFMT.h"
#include "dSFMT-src-2.2.3/dSFMT.c"

double simpi(unsigned long n){
    srand((unsigned)time(NULL));
    dsfmt_t dsfmt;
    dsfmt_init_gen_rand(&dsfmt, rand());
    unsigned long count = 0;
    double x,y;
    for(unsigned long j = 0; j < n; j++){
        x = dsfmt_genrand_close_open(&dsfmt);
        y = dsfmt_genrand_close_open(&dsfmt);
        if(x*x + y*y <= 1){
            count++;
        }
    }
    return ((double)4.0)*((double)count)/((double)n);
}
"""

filename = tempname()
filenamedl = filename * "." * Libdl.dlext

open(`gcc -Wall -O3 -march=native -xc -shared -o $filenamedl -`, "w") do f
     print(f, C_code)
end

## run(`ls -l $filenamedl`);

const simpi_gcc_dSFMT_lib = filename
simpi_gcc_dSFMT(n::Int64) = ccall((:simpi, simpi_gcc_dSFMT_lib), Float64, (Int64,), n)
simpi_gcc_dSFMT(10^5)

@time simpi_gcc_dSFMT(10^9)

### Julia

@time simpi(10^9)

f(s) = max(min(log(abs(zeta(s))), 2), -4)
x = range(0, 1, length=100)
y = range(10, 45, length=400)
s = @. x + im*y'
@time w = f.(s)
plt.figure(figsize=(8,1.5))
plt.pcolormesh(imag(s), real(s), w, cmap="gist_rainbow_r")
plt.hlines(0.5, minimum(y), maximum(y), colors="grey", lw=0.5, linestyles="dashed")
plt.xlabel("y")
plt.ylabel("x")
plt.title("log(abs(zeta(x+iy)))");

@vars a r positive=true
@vars x real=true
I = sympy.Integral(exp(-a*x^r), (x,0,oo))
sol = simplify(I.doit())
ldisp(sympy.latex(I), " = ", sympy.latex(sol))

### これは x² = ax + b を満たす x の連分数展開

function f(a, b, x, n)
    s = a + b/x
    for i in n:-1:1
        s = a+b/s
    end
    s
end

### 連分数で x² = x + 1 の解 x (黄金分割比)を数値計算

f(1,1,1,37), (1+√5)/2

### SymPyを使えば連分数を整形して表示できる.

@vars a b x
for n in 0:3
  ldisp("f(a,b,x,$n) = ", sympy.latex(f(a, b, x, n)))
end

R"""
library(ggplot2)
str(iris)
""";

R"""
p <- ggplot(iris, aes(x=Sepal.Length, y=Sepal.Width))
p + geom_point(colour="gray50", size=3) + geom_point(aes(colour=Species))
"""

n = 100
X = randn(n,2)
b = [2.0, 3.0]
y = X * b + randn(n,1)

@rput y
@rput X

R"mod <- lm(y ~ X-1)"
R"summary(mod)"

# 線分を描く函数
segment(A, B; color="black", kwargs...) = plot([A[1], B[1]], [A[2], B[2]]; color=color, kwargs...)
segment!(A, B; color="black", kwargs...) = plot!([A[1], B[1]], [A[2], B[2]]; color=color, kwargs...)
segment!(p, A, B; color="black", kwargs...) = plot!(p, [A[1], B[1]], [A[2], B[2]]; color=color, kwargs...)

# 円周を描く函数
function circle(O, r; a=0, b=2π, color="black", kwargs...)
    t = linspace(a, b, 1001)
    x(t) = O[1] + r*cos(t)
    y(t) = O[1] + r*sin(t)
    plot(x.(t), y.(t); color=color, kwargs...)
end
function circle!(O, r; a=0, b=2π, color="black", kwargs...)
    t = linspace(a, b, 1001)
    x(t) = O[1] + r*cos(t)
    y(t) = O[2] + r*sin(t)
    plot!(x.(t), y.(t); color=color, kwargs...)
end
function circle!(p, O, r; a=0, b=2π, color="black", kwargs...)
    t = linspace(a, b, 1001)
    x(t) = O[1] + r*cos(t)
    y(t) = O[1] + r*sin(t)
    plot!(p, x.(t), y.(t); color=color, kwargs...)
end

function MLIanim()
    lw = 1.0 # 太さ
    A, B, C, D = [1.5, 0], [3.5, 0], [5.5, 0], [7.5, 0]
    N = 9
    θ₀ = 2π/(2N)
    R = [
        cos(θ₀) -sin(θ₀)
        sin(θ₀)  cos(θ₀)
    ]
    V(θ) = [cos(θ), sin(θ)]
    r = 1.55*2π/(2N)

    aa = [0.5:0.025:1.5; 1.475:-0.025:0.525]
    prog = Progress(length(aa),1)
    anim = @animate for a in aa
        θ = a*2π/4
        p = plot(xlim=(-9, 9), ylim=(-9, 9))
        plot!(grid=false, legend=false, xaxis=false, yaxis=false)
        for k in 1:10
            RR = R^(2k-1)
            segment!(RR*A, RR*B,               lw=lw, color=:black)
            segment!(RR*B, RR*C,               lw=lw, color=:blue)
            segment!(RR*C, RR*D,               lw=lw, color=:red)
            segment!(RR*A, RR*(A+r*V(θ)),      lw=lw, color=:black)
            segment!(RR*A, RR*(A+r*V(-θ)),     lw=lw, color=:black)
            segment!(RR*B, RR*(B+1.5r*V(π-θ)), lw=lw, color=:black)
            segment!(RR*B, RR*(B+1.5r*V(π+θ)), lw=lw, color=:black)
            segment!(RR*C, RR*(C+2.0r*V(θ)),   lw=lw, color=:black)
            segment!(RR*C, RR*(C+2.0r*V(-θ)),  lw=lw, color=:black)
            segment!(RR*D, RR*(D+2.5r*V(π-θ)), lw=lw, color=:black)
            segment!(RR*D, RR*(D+2.5r*V(π+θ)), lw=lw, color=:black)
        end
        plot(p, size=(500, 500))
        next!(prog)
    end
    anim
end

anim = MLIanim()
gifname = "dynamic_Muller-Lyer.gif"
gif(anim, gifname, fps = 15)
sleep(0.1)
showimg("image/gif", gifname)

mixnormal(a,b) = MixtureModel([Normal(0,1), Normal(b,1)], [a, 1-a])
lpdf(a, b, y) = log(a+(1-a)*exp(b*y-b^2/2)) - y^2/2 - log(2π)/2
loglik(a, b, Y) = sum(lpdf(a, b, Y[i]) for i in 1:lastindex(Y))

function plot_lik(a₀, b₀, n; seed = 4649, bmin=-1.0, bmax=1.0)
    seed!(seed)
    Y = rand(mixnormal(a₀, b₀), n)
    L = 201
    a = range(0, 1, length=L)
    b = range(bmin, bmax, length=L)
    f(a, b) = loglik(a, b, Y)
    z = f.(a', b)
    zmax, k = findmax(z)
    #i, j = (k - 1) ÷ L + 1, mod1(k, L) # for Julia v0.6
    j, i = k.I
    z .= exp.(z .- zmax)

    sleep(0.1)
    plt.figure(figsize=(5,5))
    plt.pcolormesh(a, b, z, cmap="CMRmap")
    plt.scatter([a₀], [b₀], label="true", color="cyan")
    plt.scatter([a[i]], [b[j]], label="MLE", color="red")
    plt.legend()
    plt.xlim(0,1)
    plt.xlabel("a")
    plt.ylabel("b")
    plt.title("\$a_0 = $a₀,  b_0 = $b₀,  n = $n\$")
end

@time plot_lik(0.5, 0.1, 2^9);

@time plot_lik(0.5, 0.1, 2^12);

### Shannon情報量の収束の様子
seed!(2019)
ber = Bernoulli(1/3)
S = entropy(ber)
nmax, ntrials = 2^10, 5
n = 1:nmax
a = cumsum(rand(ber, nmax, ntrials), dims=1)./n
y = @. -(a*log(ber.p) + (1-a)*log(1-ber.p))
plot(size=(500, 300), y)
hline!([S], color=:black, ls=:dash)

seed!(181818)
X = cumsum(randn(2^14,3), dims=1)
fig = plt.figure(figsize=(6.4,4))
ax = fig.add_subplot(111, projection="3d")
ax.plot(X[:,1], X[:,2], X[:,3], lw=0.4);

function randbdr(m, n)
    k = rand(1:m+n-1)
    if k ≤ m
        return k, 1
    else
        return 1, k-m+1
    end
end

function istouching(s, i, j)
    m, n = size(s)
    s[i,j] ≠ 0 && return true
    s[ifelse(i+1>m, 1, i+1), j] ≠ 0 && return true
    s[ifelse(i-1<1, m, i-1), j] ≠ 0 && return true
    s[i, ifelse(j+1>n, 1, j+1)] ≠ 0 && return true
    s[i, ifelse(j-1<1, n, j-1)] ≠ 0 && return true
    false
end
    
function DLA_update!(s, c)
    m, n = size(s)
    i, j = randbdr(m, n)
    while !istouching(s, i, j)
        if rand() < 0.5
            i = ifelse(rand() < 0.5, ifelse(i+1>m, 1, i+1), ifelse(i-1<1, m, i-1))
        else
            j = ifelse(rand() < 0.5, ifelse(j+1>n, 1, j+1), ifelse(j-1<1, n, j-1))
        end
    end
    s[i,j] = c
    return s
end

function DLA(n, niters; seed=2019)
    seed!(seed)
    s = zeros(Int8, n, n)
    s[n÷2+1, n÷2+1] = 10
    for k in 1:niters
        c = mod(k÷(niters÷2^3), 10)+1
        DLA_update!(s, c)
    end
    s
end

function DLAseries(n, niters; nframes=200, seed=2019)
    seed!(seed)
    skip = niters ÷ nframes
    N = niters ÷ skip
    S = zeros(Int8, n, n, N+1)
    s = zeros(Int8, n, n)
    s[n÷2+1, n÷2+1] = 10
    t = 1
    S[:,:,t] = copy(s)
    for k in 1:niters
        c = mod(k÷(niters÷2^3), 10)+1
        DLA_update!(s, c)
        if mod(k, skip) == 0
            t += 1
            S[:,:,t] = copy(s)
        end
    end
    S
end

@time s = DLA(200, 4000)
@time heatmap(s, size=(400,400), color=:thermal)

@time S = DLAseries(200, 4000)
L = size(S,3)
@time anim = @animate for i in [1:2:L; L-1:-2:2]
    heatmap(S[:,:,i], size=(400,400), color=:thermal)
end

gifname = "DLA.gif"
@time gif(anim, gifname)
sleep(0.1)
showimg("image/gif", gifname)

sol2p(sol) = [sol.u[i][j] for i in 1:length(sol.u), j in 1:length(sol.u[1])÷2]
sol2q(sol) = [sol.u[i][j] for i in 1:length(sol.u), j in length(sol.u[1])÷2+1:length(sol.u[1])]
sol2tqp(sol) = (sol.t, sol2q(sol), sol2p(sol))

function H_opentoda(p, q, param=Float64[])
    sum(p.^2)/2 + sum(exp.(q[2:end] .- q[1:end-1]))
end

function solve_opentoda(q₀, p₀; integrator=Yoshida6(), Δt=0.01, t₁=10.0)
    prob = HamiltonianProblem(H_opentoda, p₀, q₀, (0.0, t₁))
    solve(prob, integrator, dt=Δt)
end

function solve_discopentoda(X₀; t₁=10)
    solX = Array{typeof(X₀),1}(undef, t₁+1)
    solX[1] = X₀
    for k in 2:t₁+1
        C = cholesky(solX[k-1])
        solX[k] = C.U*C.L
    end
    return solX
end

function Bmat(L)
    B = -one(L)
    B[1:end-1, 2:end] -= @view B[1:end-1, 1:end-1]
    B[end,:] = ones(size(B)[2])
    return B
end

function L2q(L)
    n = size(L)[1]
    qq = zeros(eltype(L), n)
    for i in 1:n-1
        qq[i] = 2log(L[i,i+1])
    end
    q = Bmat(L)\qq
end
L2qp(L) = (L2q(L), diag(L))

function qp2L(q, p)
    L = Matrix(Diagonal(p))
    for i in 1:length(q)-1
        L[i,i+1] = L[i+1,i] = exp(0.5*(q[i+1]-q[i]))
    end
    return L
end

X2qp(X) = L2qp(log(X))
qp2X(q,p) = exp(qp2L(q,p))

solX2q(solX) = [X2qp(solX[i])[1][j] for i in eachindex(solX), j in 1:size(solX[1])[1]]
solX2p(solX) = [X2qp(solX[i])[2][j] for i in eachindex(solX), j in 1:size(solX[1])[1]]
solX2qp(solX) = (solX2q(solX), solX2p(solX));

colorlist = [
    "#1f77b4", "#ff7f0e", "#2ca02c", "#d62728", "#9467bd",
    "#8c564b", "#e377c2", "#7f7f7f", "#bcbd22", "#17becf",
]
cc(k) = colorlist[mod1(k, length(colorlist))]

function plot_opentoda(q₀, p₀; t₁=10)
    q₀₀ = q₀ .- mean(q₀)
    d = length(p₀)
    
    sol = solve_opentoda(q₀₀, p₀, t₁=Float64(t₁))
    t, q, p = sol2tqp(sol)
    
    X₀ = qp2X(q₀₀, p₀)
    solX = solve_discopentoda(X₀, t₁=t₁)
    qX, pX = solX2qp(solX)
    
    P1 = plot(legend=false)
    for j in 1:d
        plot!(t, q[:,j], color=cc(j), lw=1)
        scatter!(0:t₁, qX[:,j],  color=cc(j), markersize=2)
    end
    xlabel!("t")
    ylabel!("q")
    title!("discrete and continuous open Toda: q_i", titlefontsize=10)

    P2 = plot(legend=false)
    for j in 1:d
        plot!(t, p[:,j], color=cc(j), lw=1)
        scatter!(0:t₁, pX[:,j],  color=cc(j), markersize=2)
    end
    xlabel!("t")
    ylabel!("p")
    title!("discrete and continuous open Toda: p_i", titlefontsize=10)
    
    plot(P1, P2, size=(700, 300))
end

q₀ = Float64[4, 3, 2, 1]
p₀ = Float64[-2, -1, 1, 2]
plot_opentoda(q₀, p₀)

q₀ = Float64[5, 2, 4, 1]
p₀ = Float64[-3, -2, 1, 4]
plot_opentoda(q₀, p₀)

q₀ = Float64[10, 2, 6, 1]
p₀ = Float64[-3, -2, 1, 4]
plot_opentoda(q₀, p₀)

n = 2^10
seed!(2018)
X = randn(n,n)/√n # すべての成分が平均0, 分散1/nの正規分布に従う
@time λ, U = eigen(X)

plt.figure(figsize=(4, 4))
plt.scatter(real(λ), imag(λ), s=10.0)
plt.grid(ls=":")
plt.title("Circular law");

n = 2^10
seed!(2018)
X = (2rand(n,n) .- 1) * √(3/n) # [-√(3/n), √(3/n)]上の一様分布. 分散は1/nになる.
@time λ, U = eigen(X)

plt.figure(figsize=(4, 4))
plt.scatter(real(λ), imag(λ), s=10.0)
plt.grid(ls=":")
plt.title("Circular law");

n = 2^10
seed!(2019)
X = Symmetric(randn(n,n)) # 標準正規分布の場合
@time λ, U = eigen(X)

a = λ/√n
x = range(-2, 2, length=200)
f(x) = 1/(2π)*sqrt(4-x^2)

plt.figure(figsize=(5, 2.7))
plt.hist(a, normed=true, bins=50, alpha=0.5, label="normaized eigen values")
plt.plot(x, f.(x), label="\$y = (2/π)\\sqrt{1-x^2}\$", color="red", ls="--")
plt.grid(ls=":")
plt.legend()
plt.title("Semicircle law");

n = 2^10
seed!(2019)
X = Symmetric((2*√3)*rand(n,n) .- √3) # 分散1の一様分布の場合
@time λ, U = eigen(X)

a = λ/√n
x = range(-2, 2, length=200)
f(x) = 1/(2π)*sqrt(4-x^2)

plt.figure(figsize=(5, 2.7))
plt.hist(a, normed=true, bins=50, alpha=0.5, label="normaized eigen values")
plt.plot(x, f.(x), label="\$y = (2/π)\\sqrt{1-x^2}\$", color="red", ls="--")
plt.grid(ls=":")
plt.legend()
plt.title("Semicircle law");

@vars x
I = sympy.Integral(√(1-x^2), (x,-1,1))
ldisp(sympy.latex(I), " = ", sympy.latex(I.doit()))

@vars θ
J = sympy.Integral(sin(θ)^2, (θ, 0, PI))
ldisp(sympy.latex(J), " = ", sympy.latex(J.doit()))

function nrationalpoints_naive(f, p)
    @distributed (+) for x in 0:p-1
        s = 0
        for y in 0:p-1
            s += ifelse(mod(f(x,y),p) == 0, 1, 0)
        end
        s
    end
end

function plot_SatoTate_naive(f; figtitle="Sato-Tate conjecture", N=2^12)
    P = primes(N)
    @show N, length(P)
    @time S = nrationalpoints_naive.(f, P) .+ 1 # "+1" は無限遠点の個数
    plot_SatoTate(P, S; figtitle=figtitle)
end

function nrationalpoints_legendre(g, p)
    @distributed (+) for x in 0:p-1
        l = legendresymbol(mod(g(x),p), p)
        ifelse(l == 1, 2, ifelse(l == -1, 0, 1))
    end
end

function plot_SatoTate_legendre(f; figtitle="Sato-Tate conjecture", N=2^12)
    P = primes(N)
    @show N, length(P)
    @time S = nrationalpoints_legendre.(f, P) .+ 1 # "+1" は無限遠点の個数
    plot_SatoTate(P, S; figtitle=figtitle)
end

function plot_SatoTate(P, S; figtitle="Sato-Tate conjecture")
    a = (P .+ 1) - S

    @show count(abs.(a) .> 2sqrt.(P)) # Weil予想の確認
    X = a ./ (2sqrt.(P)) # -1 から 1 の区間に入るように正規化
    θ = acos.(X)

    x = range(-1, 1, length=200)
    g(x) = (2/π)*sqrt(1-x^2) # 半円則

    t = range(0, π, length=200)
    h(t) = (2/π)*sin(t)^2 # sin^2 分布

    sleep(0.1)
    plt.figure(figsize=(8,3))
    
    plt.subplot(121)
    plt.hist(X, normed=true, bins=21, alpha=0.5, label="\$a_p/(2\\sqrt{p})\$")
    plt.plot(x, g.(x), color="red", ls="--", label="\$y=(2/\\pi)\\sqrt{1-x^2}\$")
    plt.xlabel("\$x\$")
    plt.grid(ls=":")
    plt.legend(fontsize=9)
    plt.title(figtitle, fontsize=10)
    
    plt.subplot(122)
    plt.hist(θ, normed=true, bins=21, alpha=0.5, label="\$\\arccos(a_p/(2\\sqrt{p}))")
    plt.plot(t, h.(t), color="red", ls="--", label="\$y=(2/\\pi)\\sin^2\\theta\$")
    plt.xlabel("\$\\theta\$")
    plt.grid(ls=":")
    plt.legend(fontsize=8)
    plt.title(figtitle, fontsize=10)

    plt.tight_layout()
end

addedprocs = addprocs(4)
@everywhere f(x,y) = y^2 + y - (x^3-x^2)
N = 2^12
figtitle = "Sat-Tate conj. for  \$y^2 + y = x^3-x^2\$,  \$p < 2^{12}\$"
plot_SatoTate_naive(f, figtitle=figtitle, N=N)
rmprocs(addedprocs);

addedprocs = addprocs(4)
@everywhere using Combinatorics: legendresymbol
@everywhere g(x) = x^3+x+1
N = 2^14
figtitle = "Sat-Tate conj. for  \$y^2 = x^3+x+1\$,  \$p < 2^{14}\$"
plot_SatoTate_legendre(g, figtitle=figtitle, N=N)
rmprocs(addedprocs);

using FileIO: load, save
using LibSndFile: LibSndFile

## downloaded from https://freewavesamples.com/casio-mt-600-harpsichord-c3
samplesound1 = load("sounds/Casio-MT-600-Harpsichord-C3.wav")

## downloaded from http://music.futta.net/mp3.html
samplesound2 = load("sounds/futta-dream.wav")

# plotlyjs()
# x = range(-3, 3, length=300)
# y = range(-3, 3, length=300)
# z = @. exp(-(x'^2+y^2))
# clibrary(:misc)
# P = surface(x, y, z, colorbar=false, size=(600, 500), color=:rainbow)
# display(P)
# gr();