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)