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", """""") 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 #include 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 #include #include #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();