using LinearAlgebra, Plots
bucket = "https://raw.githubusercontent.com/nmfsc/data/master/";

packages = ["LinearAlgebra","Plots","Images","DelimitedFiles","SparseArrays",
  "Graphs","GraphPlot","FFTW","Primes","BenchmarkTools","DifferentialEquations",
  "Optim","Dierckx","Interpolations","SpecialFunctions","Roots","OffsetArrays",
  "Wavelets","Interact","LsqFit","GLM","DataFrames","Flux","FastGaussQuadrature",
  "Sundials","SciMLOperators","CSV","ColorSchemes","JuMP","GLPK","MLDatasets",
  "Arpack","ProgressMeter","IterativeSolvers","Polynomials","ParameterizedFunctions",
  "ImageShow","NLsolve","DiffEqDevTools","LinRegOutliers","ModelingToolkit","NonlinearSolve"];
#import Pkg; Pkg.add(packages)

hilbert(n)  = [1/(i+j-1) for i=1:n, j=1:n]

using Images
[Gray.(1 .- abs.(hilbert(n)\hilbert(n))) for n ∈ (10,15,20,25,50)]

function gaussian_elimination(A,b)
  n = size(A,1)
  for j in 1:n
    A[j+1:n,j] /= A[j,j]
    A[j+1:n,j+1:n] -= A[j+1:n,j:j].*A[j:j,j+1:n]
  end
  for i in 2:n
    b[i:i] -= A[i:i,1:i-1]*b[1:i-1]
  end
  for i in n:-1:1
    b[i:i] = ( b[i] .- A[i:i,i+1:n]*b[i+1:n] )/A[i,i]
  end
  return b
end

A = rand(8,8); b = rand(8,1)
[A\b gaussian_elimination(A,b)]

ϵ = 1e-15; A = [1 1 2;2 2+ϵ 0; 4 14 4]; b = [-5;10;0.0]
A\b, gaussian_elimination(A,b)

function row_reduce!(tableau)
  (i,j) = get_pivot(tableau)
  G = tableau[i:i,:]/tableau[i,j]
  tableau[:,:] -= tableau[:,j:j]*G
  tableau[i,:] = G
end

function get_pivot(tableau)
  j = argmax(tableau[end,1:end-1])
  a, b = tableau[1:end-1,j], tableau[1:end-1,end]
  k = findall(a.>0)
  i = k[argmin(b[k]./a[k])]
  return(i,j) 
end

function simplex(c,A,b)
  (m,n) = size(A)
  tableau = [[A I b] ; [c' zeros(1,m) 0]]
  while (any(tableau[end,1:n].>0))
    row_reduce!(tableau)
  end
  p = findall(tableau[end,1:n].==0)
  x = zeros(n,1)
  [x[i]=tableau[:,i]⋅tableau[:,end]  for i∈p]
  z = -tableau[end,end]
  y = -tableau[end,n.+(1:m)]
  return((z=z, x=x, y=y))
end

A = [2 0 1;4 1 2;1 1 0];  c = [2;1;1];  b = [3;2;1]
solution = simplex(c,A,b)

using DelimitedFiles, SparseArrays
function get_adjacency_matrix(filename)
  ij = readdlm(download(bucket*filename*".csv"),',',Int)
  sparse(ij[:,1],ij[:,2],one.(ij[:,1]))
end

using Graphs, GraphPlot
g = get_adjacency_matrix("dolphins") |> Graph
gplot(g,layout=spectral_layout)

gplot(g,layout=spring_layout)

gplot(g,layout=circular_layout)

using SparseArrays
function revised_simplex(c,A,b)
  (m,n) = size(A)
  ξ = i -> (z=spzeros(m);z[i]=1;z)
  N = Vector(1:n); B = Vector(n .+ (1:m))
  A = [A sparse(I, m, m)]
  ABinv = sparse(I, m, m)
  c = [c;spzeros(m,1)]
  while(true)
    j = findfirst(x->x>0,(c[N]'.-(c[B]'*ABinv)*A[:,N])[:])
    if isnothing(j); break; end
    q = ABinv*A[:,N[j]]
    k = findall(q.>0)
    i = k[argmin(ABinv[k,:]*b./q[k])]
    B[i], N[j] = N[j], B[i]
    ABinv -=  ((q - ξ(i))/q[i])*ABinv[i:i,:]
  end
  i = findall(B.≤n)
  x = zeros(n,1)
  x[B[i]] = ABinv[i,:]*b
  return((z=c[1:n]'*x, x=x, y=c[B]'*ABinv))
end

A = [2 0 1;4 1 2;1 1 0];  c = [2;1;1];  b = [3;2;1]
solution = simplex(c,A,b)

using DelimitedFiles
T = readdlm(download(bucket*"sherlock.csv"), '\t')[:,2]
n = length(T)
A = [ones(n,1) log.(1:n)]
B = log.(T)
c = A\B
print("ordinary least squares:\n"*string(c)*"\n")

zipf(x,c) = exp(c[1])*x.^c[2]
scatter(T,xaxis=:log,yaxis=:log,ma=0.4,msw=0)
plot!([1,n],zipf([1,n],c),w=3)

function constrained_lstsq(A,b,C,d)
  x = [A'*A C'; C zeros(size(C,1),size(C,1))]\[A'*b;d]
  x[1:size(A,2)]
end

function tls(A,B)
  n = size(A,2)
  V = svd([A B]).V
  -V[1:n,n+1:end]/V[n+1:end,n+1:end]
end

A = [2 4; 1 -1; 3 1; 4 -8]; b = [1; 1; 4; 1];
xₒₗₛ = A\b; xₜₗₛ = tls(A,b)

print("ordinary least squares:"*string(xₒₗₛ))
print("\ntotal least squares:"*string(xₜₗₛ))

using Images
A = load(download(bucket*"laura.png")) .|> Gray
U, σ, V = svd(A);

k = 20
Aₖ = U[:,1:k] * Diagonal(σ[1:k]) * V[:,1:k]' .|> Gray
[A Aₖ]

norm(A-Aₖ) ≈ norm(σ[k+1:end])

ϵ² = 1 .- cumsum(σ)/sum(σ); scatter(ϵ²,xaxis=:log)

function nmf(X,p=6)
  W = rand(Float64, (size(X,1), p))
  H = rand(Float64, (p,size(X,2)))
  for i in 1:50
    W = W.*(X*H')./(W*(H*H') .+ (W.≈0))
    H = H.*(W'*X)./((W'*W)*H .+ (H.≈0))
  end
  (W,H)
end

using Images
A = Gray.(load(download(bucket*"laura.png")))
W,H = nmf(Float64.(A),20)
[A Gray.(W*H)]

function condeig(A)
  Sᵣ = eigvecs(A)
  Sₗ = inv(Sᵣ')
  Sₗ ./= sqrt.(sum(abs.(Sₗ.^2), dims=1))
  1 ./ abs.(sum(Sᵣ.*Sₗ, dims=1))
end

condeig(rand(4,4))

H = [0 0 0 0 1; 1 0 0 0 0; 1 0 0 0 1; 1 0 1 0 0; 0 0 1 1 0]
v = all(H.==0,dims=1)
H = H ./ (sum(H,dims=1)+v) 
n = size(H,1) 
d = 0.85
x = ones(n,1)/n
for i in 1:9
  x = d*(H*x) .+ d/n*(v*x)  .+ (1-d)/n
end 
x

using  FFTW, Images, ColorSchemes
perfectshuffle= A -> [A[1:2:end,:];A[2:2:end,:]]
F = real.(fft(I(128),1)); PF = F
for i = 1:6
  F = perfectshuffle(F); PF = [PF F]
end
get(colorschemes[:gray1], PF, :extrema)

function fftx2(c)
  n = length(c)
  ω = exp(-2im*π/n)
  if mod(n,2) == 0 
    k = collect(0:n/2-1)
    u = fftx2(c[1:2:n-1])
    v = (ω.^k).*fftx2(c[2:2:n]) 
    return([u+v; u-v])
  else 
    k = collect(0:n-1)
    F = ω.^(k*k')
    return(F*c)
  end 
end

ifftx2(y) = conj(fftx2(conj(y)))/length(y);

function fasttoeplitz(c,r,x)
  n = length(x)
  Δ = nextpow(2,n) - n
  x₁ = [c; zeros(Δ); r[end:-1:2]]
  x₂ = [x; zeros(Δ+n-1)]
  ifftx2(fftx2(x₁).*fftx2(x₂))[1:n]
end

using FFTW, Primes, Plots
N = 10000
smooth(n,N) = (1:N)[all.(x->x<=n,factor.(Set,1:N))]
t₁ = [(x = randn(n); (n,@elapsed(fft(x)))) for n∈primes(N)]
t₂ = [(x = randn(n); (n,@elapsed(fft(x)))) for n∈smooth(5,N)]
plot(t₁,label="prime"); plot!(t₂,label="5-smooth")
plot!(ylims=(0,t₁[end][2]*3))

function bluestein_fft(x)
  n = length(x)
  ω = exp.((1im*π/n)*(0:n-1).^2)
  ω.*fasttoeplitz(conj(ω),conj(ω),ω.*x)
end

using Primes
function primitiveroot(n)
  ϕ = n - 1
  p = factor(Set,ϕ)
  for r = 2:n-1
    all( [powermod(r, ϕ÷pᵢ, n) for pᵢ∈p] .!= 1 ) && return(r)
  end
end

function rader_fft(x)
  n = length(x)
  r = primitiveroot(n)
  P₊ = powermod.(r, 0:n-2, n)
  P₋ = circshift(reverse(P₊),1)
  ω = exp.((2im*π/n)*P₋)
  c = x[1] .+ ifft(fft(ω).*fft(x[2:end][P₊]))
  [sum(x); c[reverse(invperm(P₋))]]    
end

using FFTW
dst(x) = FFTW.r2r(x,FFTW.RODFT00)
idst(x) = dst(x)/(2^ndims(x)*prod(size(x).+1))

n = 50; x = (1:n)/(n+1); Δx = 1/(n+1) 
v = 2 .- 2cos.(x*π) 
λ = [v[i]+v[j]+v[k] for i∈1:n, j∈1:n, k∈1:n]./Δx^2
f = [(x-x^2)*(y-y^2) + (x-x^2)*(z-z^2)+(y-y^2)*(z-z^2) 
    for x∈x,y∈x,z∈x]
u = idst(dst(f)./λ);

norm(u - [(x-x^2)*(y-y^2)*(z-z^2)/2 for x∈x,y∈x,z∈x])

using Images, FFTW, SparseArrays
function dctcompress(A,d)
  B = dct(Float64.(A))
  idx = d*prod(size(A)) |> floor |> Int
  tol = sort!(abs.(B[:]),rev=true)[idx]
  B₀ = droptol!(sparse(B),tol)
  idct(B₀) |> clamp01! .|> Gray
end

A = Gray.(load(download(bucket*"laura_square.png")))
[A dctcompress(A,0.05)]

bitstring(Float64(π))

function Q_rsqrt(x::Float32)
  i = reinterpret(Int32,Float32(x)) 
  i = Int32(0x5f3759df) - (i>>1)   
  y = reinterpret(Float32,i)   
  y = y * (1.5f0 - (0.5f0 * x * y * y))     
end

Q_rsqrt(Float32(0.01))

using BenchmarkTools
@btime Q_rsqrt(Float32(x)) setup=(x=rand()) seconds=3
@btime 1/sqrt(x) setup=(x=rand()) seconds=3

a = 77617; b = 33096
333.75*b^6+a^2*(11*a^2*b^2-b^6-121*b^4-2)+5.5*b^8+a/(2*b)

function bisection(f,a,b,tolerance)
  while abs(b-a) > tolerance
    c = (a+b)/2
    sign(f(c)) == sign(f(a)) ? a = c : b = c
  end
  (a+b)/2
end

bisection(x->sin(x),2,4,1e-14)

function escape(n, c=0, z=0)
  for k in 1:n
    z = z^2 + c
    abs2(z) > 4 && return k
  end
  return n
end

function mandelbrot(bb, xn=800, n=200, z=0)
  yn = (xn*(bb[4]-bb[2]))÷(bb[3]-bb[1]) |> Int
  c = LinRange(bb[1],bb[3],xn)' .+ im*LinRange(bb[4],bb[2],yn)
  [escape(n,c,z) for c∈c]
end

using Images
M = mandelbrot([-0.172, 1.0228, -0.1494, 1.0443])
save("mandelbrot.png",1 .- M./maximum(M))
load("mandelbrot.png")

function julia(bb, xn=800, n=200, c=-0.123+0.745im)
  yn = (xn*(bb[4]-bb[2]))÷(bb[3]-bb[1]) |> Int
  z = LinRange(bb[1],bb[3],xn)' .+ im*LinRange(bb[4],bb[2],yn)
  [escape(n,c,z) for z∈z]
end

J = julia([-2, -2, 2, 2])
save("julia.png",1 .- J./maximum(J))
load("julia.png")

using DifferentialEquations
f = z -> ((x,y)=tuple(z...); [x^3-3x*y^2-1; y^3-3x^2*y] )
df = (z,p,_) -> ((x,y)=tuple(z...);
  -[3x^2-3y^2 -6x*y; -6x*y 3y^2-3x^2]\p)
z₀ = [1,1]
sol = solve(ODEProblem(df,z₀,(0,1),f(z₀)))
sol.u[end]

function plot_trace(x,f)
  contour(-2:0.02:1, -0.25:0.02:3.25,(x,y) -> cbrt(f([x,y])), color=:red, colorbar=:none)
  plot!(x[1,:],x[2,:],color=:black, linewidth=2, marker=:o, legend=:none)
end

∇f = x-> [-2(1-x[1])-400x[1]*(x[2]-x[1]^2); 200(x[2]-x[1]^2)]
x = [-1.8,3.0]; p = [0,0]; α =  0.001; β = 0.9
s = x
for i = 1:100
  p = -∇f(x) + β*p
  x += α*p
  append!(s,x)
end
f = x -> (1-x[1])^2 + 100(x[2] - x[1]^2)^2
plot_trace(reshape(s, 2, :),f)

using Optim
x₀ = [-1.8,3.0]
options = Optim.Options(store_trace = true,extended_trace = true)
result = optimize(f, x₀, BFGS(),options)
s = Optim.x_trace(result); s = vcat(s...)
plot_trace(reshape(s, 2, :),f)

∇f = x -> [-2(1-x[1])-400x[1]*(x[2]-x[1]^2); 200(x[2]-x[1]^2)]
x = [-1.8,3.0]; α =  0.001; p = [0,0]; β = 0.9
for i = 1:100
  p = -∇f(x) + β*p
  x += α*p
end

using Optim
f = x ->  (1-x[1])^2 + 100(x[2] - x[1]^2)^2
x₀ = [-1.8,3.0]
result = optimize(f, x₀, GradientDescent())
x = Optim.minimizer(result)

function spline_natural(xᵢ,yᵢ)
  h = diff(xᵢ)
  γ = 6*diff(diff(yᵢ)./h)
  α = h[2:end-1]
  β = 2(h[1:end-1]+h[2:end])
  [0;SymTridiagonal(β,α)\γ;0]
end

function evaluate_spline(xᵢ,yᵢ,m,n)
  h = diff(xᵢ)
  B = yᵢ[1:end-1] .- m[1:end-1].*h.^2/6
  A = diff(yᵢ)./h - h./6 .*diff(m)
  x = range(minimum(xᵢ),maximum(xᵢ),length=n+1)
  i = sum(xᵢ'.≤x,dims=2)
  i[end] = length(xᵢ)-1 
  y = @. (m[i]*(xᵢ[i+1]-x)^3 + m[i+1]*(x-xᵢ[i])^3)/(6h[i]) +
     A[i]*(x-xᵢ[i]) + B[i]
  return(x,y)
end

x = LinRange(0,1,8); y = rand(8)
m = spline_natural(x,y)
X,Y = evaluate_spline(x,y,m,100)
scatter(x,y); plot!(X,Y)

g = x-> @. max(1-abs(3-x),0)
xᵢ = 0:5; yᵢ = g(xᵢ)
x = LinRange(0,5,101);

using Dierckx, Plots
spline = Spline1D(xᵢ,yᵢ)
plot(x,spline(x)); plot!(x,g(x)); scatter!(xᵢ,yᵢ)

using Interpolations
method = BSpline(Cubic(Natural(OnGrid())))
spline =  scale(interpolate(yᵢ, method), xᵢ)
plot(x,spline(x)); plot!(x,g(x)); scatter!(xᵢ,yᵢ)

bernstein(n,t) = @. binomial(n,0:n)'*t^(0:n)'*(1-t)^(n:-1:0)'

n = 3
t = LinRange(0,1,100)
p = rand(n+1,2)
z = bernstein(n,t)*p
plot(p[:,1],p[:,2],marker=:o,opacity=0.3); 
plot!(z[:,1],z[:,2],legend=:none)

n = 20; f = t -> sqrt(t)
y = bernstein(n,t)*f.(LinRange(0,1,n+1))
plot(t,y); plot!(t,f.(t))

function legendre(x,n)
  n==0 && return(one.(x))
  n==1 && return(x)
  x.*legendre(x,n-1) .- (n-1)^2/(4(n-1)^2-1)*legendre(x,n-2)
end

x = LinRange(-1,1,100)
plot()
for n in 0:4
    plot!(x,legendre(x,n)) 
end
current()

function chebdiff(n)
  x = -cos.(LinRange(0,π,n))
  c = [2;ones(n-2);2].*(-1).^(0:n-1)
  D = c./c'./(x .- x' + I)
  D - Diagonal([sum(D,dims=2)...]), x  
end

n = 15
D,x = chebdiff(n)
u = exp.(-4x.^2)
plot(x,D*u,marker=:o)
t = LinRange(-1,1,200)
plot!(t,-8t.*exp.(-4t.^2))

B = zeros(1,n); B[1] = 1
plot(x,[B;D]\[2;u],marker=:o)

n = 15; k² = 256
D,x = chebdiff(n)
L = (D^2 - k²*Diagonal(x))
L[[1,n],:] .= 0; L[1,1] = L[n,n] = 1 
y = L\[2;zeros(n-2);1]
plot(x,y,marker=:o)

using SpecialFunctions
a = [airyai(-∛k²) airybi(-∛k²);airyai(∛k²) airybi(∛k²)]\[2;1]
airy(x) = a[1]*airyai(∛k²*x) + a[2]*airybi(∛k²*x)
scatter(x,y)
plot!(t,airy.(t))

λ² = 16^2; N = 6:2:60; ϵ = []
for n ∈ N 
  D,x = chebdiff(n)
  L = (D^2 - k²*Diagonal(x))
  L[[1,n],:] .= 0; L[1,1] = L[n,n] = 1 
  y = L\[2;zeros(n-2);1]
  append!(ϵ,norm(y - airy.(x),Inf))
end
plot(N,ϵ,yaxis=:log,marker=:o,legend=:none)

using OffsetArrays
function scaling(c,ϕₖ,n)
  m = length(c)-1; ℓ = 2^n
  ϕ = OffsetVector(zeros(3*m*ℓ),-m*ℓ) 
  k = (0:m)*ℓ
  ϕ[k] = ϕₖ
  for j = 1:n
    for i = 1:m*2^(j-1)
      x = (2i-1)*2^(n-j)            
      ϕ[x] = c ⋅ ϕ[2x .- k]
    end
  end
  ((0:m*ℓ-1)/ℓ, ϕ[0:m*ℓ-1])
end

c = [1+√3, 3+√3, 3-√3, 1-√3]/4
z = [0, 1+√3, 1-√3, 0]/2
(x,ϕ) = scaling(c, z, 8)
plot(x,ϕ)

ψ = zero(ϕ); n = length(c)-1; ℓ = length(ψ)÷2n
for k∈0:n
  ψ[(k*ℓ+1):(k+n)*ℓ] += (-1)^k*c[n-k+1]*ϕ[1:2:end]
end
plot(x,ψ)

using Wavelets, Images
img = float.(Gray.(load(download(bucket*"laura_square.png"))))
B = dwt(Float32.(img), wavelet(WT.haar))
[img Gray.(1 .- clamp01.(sqrt.(abs.(B))))]

function dwtfilter(channels,wt,k)  
  map(1:3) do i
    A = dwt(channels[i,:,:], wavelet(wt))
    threshold!(A, HardTH(), k) 
    clamp01!(idwt(A, wavelet(wt)))
  end
end

using Interact, Wavelets, Images
img = float.(load(download(bucket*"laura_square.png")))
channels = channelview(img)
func = Dict("Haar"=>WT.haar, "D₄"=>WT.db4, "Coiflet"=>WT.coif4)
@manipulate  for wt in togglebuttons(func; label="Transform"),
  k in slider(0:0.05:1.0; value = 0, label = "Threshold")
  colorview(RGB,dwtfilter(channels,wt,k)...)
end

function jacobian(f,x,c)
  J = zeros(length(x),length(c))
  for k in (n = 1:length(c))
    J[:,k] .= imag(f(x,c+1e-8im*(k .== n)))/1e-8
  end
  return J
end

function gauss_newton(f,x,y,c)
  r = y - f(x,c) 
  for j = 1:100
    G = jacobian(f,x,c)
    M = G'*G
    c += (M+Diagonal(M))\(G'*r)
    norm(r-(r=y-f(x,c))) < 1e-12 && return(c)
  end
  print("Gauss-Newton did not converge.")
end

f = (x,c) -> @. c[1]*exp(-c[2]*(x-c[3])^2) +
         c[4]*exp(-c[5]*(x-c[6])^2)
x = 8*rand(100)
y = f(x,[1, 3, 3, 2, 3, 6]) + 0.1*randn(100);

c₀ = [2, 0.3, 2, 1, 0.3, 7]
c = gauss_newton(f,x,y,c₀)

if !isnothing(c)
  X = LinRange(0,8,200)
  scatter(x,y,opacity=0.5)
  plot!(X,f(X,c))
end

using LsqFit
cf = curve_fit(f,x,y,c₀)
c, r = cf.param, cf.resid
scatter(x,y,opacity=0.5)
plot!(X,f(X,c))

σ = x -> @. 1/(1+exp(-x))
x = rand(30); y = ( rand(30) .< σ(10x.-5) );

X, w = [one.(x) x], zeros(2,1)
for i=1:10
  S = σ(X*w).*(1 .- σ(X*w))
  w += (X'*(S.*X))\(X'*(y - σ(X*w)))
end

scatter(x,y)
t = LinRange(minimum(x),maximum(x),50)
plot!(t,σ([one.(t) t]*w))

using GLM, DataFrames
data = DataFrame(X=x, Y=y)
logit = glm(@formula(Y ~ X), data, Bernoulli(), LogitLink())

N = 100; θ = LinRange(0,π,N)'
x = cos.(θ); x̃ = [one.(x);x]
y = sin.(θ) + 0.05*randn(1,N)

n = 20; W₁ = rand(n,2); W₂ = randn(1,n)
ϕ = x -> max.(x,0)
dϕ = x -> (x.>0)

α = 0.1
for epoch = 1:1000
  ŷ = W₂ * ϕ(W₁*x̃)
  ∂L∂y = 2(ŷ-y)/N
  ∂L∂W₁ = dϕ(W₁*x̃) .* (W₂' * ∂L∂y) * x̃'
  ∂L∂W₂ = ∂L∂y * ϕ(W₁*x̃)'
  W₁ -= 0.1α * ∂L∂W₁ 
  W₂ -= α * ∂L∂W₂ 
end

scatter(x̃[2,:],y',opacity=0.3)
x̂ = LinRange(-1.2,1.2,200)'; x̂ = [one.(x̂);x̂]; ŷ = W₂ * ϕ(W₁*x̂)
plot!(x̂[2,:],ŷ',width=3)

n = 20; W₁ = rand(n,2); W₂ = randn(1,n)
ϕ = x ->  @. 1 / (1 + exp(-x))
dϕ = x -> @. ϕ(x)*(1-ϕ(x))
α = 0.1
for epoch = 1:10000
  ŷ = W₂ * ϕ(W₁*x̃)
  L = norm(ŷ - y)
  ∂L∂y = (ŷ-y)/L
  ∂L∂W₁ = dϕ(W₁*x̃) .* (W₂' * ∂L∂y) * x̃'
  ∂L∂W₂ = ∂L∂y * ϕ(W₁*x̃)'
  W₁ -= 0.1α * ∂L∂W₁ 
  W₂ -= α * ∂L∂W₂ 
end
scatter(x̃[2,:],y',opacity=0.3)
x̂ = LinRange(-1.2,1.2,200)'; x̂ = [one.(x̂);x̂]; ŷ = W₂ * ϕ(W₁*x̂)
plot!(x̂[2,:],ŷ',width=3)

using Flux
model = Chain(Dense(1, n, relu), Dense(n, 1))
loss(x,y) = Flux.Losses.mse(model(x), y)
parameters = Flux.params(model)
data = [(x,y)]
optimizer = Descent(0.1)
for epochs = 1:10000
  Flux.train!(loss, parameters, data, optimizer)
end

scatter(x',y',ma=0.4,msw=0)
plot!(x',model(x)',width=3)

δ = [-1,0,1,2]
n = length(δ)
V = δ.^(0:n-1)' .// factorial.(0:n-1)'
C = inv(V)

function richardson(f,x,m,n=m)
  n > 0 ? 
  (4^n*richardson(f,x,m,n-1) - richardson(f,x,m-1,n-1))/(4^n-1) :
  ϕ(f,x,2^m)  
end

ϕ = (f,x,n) -> (f(x+1/n) - f(x-1/n))/(2/n)

richardson(x->sin(x),0,4)

struct Dual
  value
  deriv
end

Dual(x) = Dual(x,1)
value(x) = x isa Dual ? x.value : x
deriv(x) = x isa Dual ? x.deriv : 0

Base.:+(u, v) = Dual(value(u)+value(v), deriv(u)+deriv(v))
Base.:-(u, v) = Dual(value(u)-value(v), deriv(u)-deriv(v))
Base.:*(u, v) = Dual(value(u)*value(v), 
   value(u)*deriv(v)+value(v)*deriv(u))
Base.:sin(u)  =  Dual(sin(value(u)), cos(value(u))*deriv(u))

x = Dual(0)
y = x + sin(x)

x₁ = Dual(2,[1 0])
x₂ = Dual(π,[0 1])
y₁ = x₁*x₂ + sin(x₂)
y₂ = x₁*x₂ - sin(x₂)
value(y₁),value(y₁),deriv(y₁),deriv(y₂)

function trapezoidal(f,x,n) 
  F = f.(LinRange(x[1],x[2],n+1))
  (F[1]/2 + sum(F[2:end-1]) + F[end]/2)*(x[2]-x[1])/n
end

ϕ = (f,x,n) -> trapezoidal(f,x,n)

richardson(x->sin(x),[0,π/2],4)

n = [floor(Int,10^y) for y in LinRange(1, 2, 10)]
error = zeros(10,7)
f = (x,p) -> x + x.^p.*(2-x).^p
for p ∈ 1:7,
  S = trapezoidal(x->f(x,p),[0,2],10^6)
  for i ∈ 1:length(n)
    Sₙ = trapezoidal(x->f(x,p),[0,2],n[i])
    error[i,p] =  abs(Sₙ - S)/S
  end
end
slope = ([log.(n) one.(n)]\log.(error))[1,:]
info = .*(string.((1:7)'),": slope=",string.(round.(slope')))
plot(n,error, xaxis=:log, yaxis=:log, labels = info)

using FFTW, LinearAlgebra
function clenshaw_curtis(f,n)
  x = cos.(π*(0:n)'/n)
  w = zeros(n+1,1); w[1:2:n+1] = 2 ./ (1 .- (0:2:n).^2)
  1/n * dctI(f.(x)) ⋅ w
end

function dctI(f) 
  g = FFTW.r2r!([f...],FFTW.REDFT00)
  [g[1]/2; g[2:end-1]; g[end]/2]  
end

clenshaw_curtis(x -> cos(8x^2),20)

function dctI(f)
  n = length(f)
  g = real(fft(f[[1:n; n-1:-1:2]]))
  [g[1]/2; g[2:n-1]; g[n]/2]  
end

function gauss_legendre(n)
  a = zeros(n)  
  b = (1:n-1).^2 ./ (4*(1:n-1).^2 .- 1)
  𝟙² = 2
  λ, v = eigen(SymTridiagonal(a, sqrt.(b))) 
  (λ, 𝟙²*v[1,:].^2)
end

n = 8
f = x -> cos(x)*exp(-x^2);
nodes, weights = gauss_legendre(n)
f.(nodes) ⋅ weights

using FastGaussQuadrature
nodes, weights = gausslegendre(n)
f.(nodes) ⋅ weights

r = exp.(2im*π*(0:.01:1))
plot(@. (1.5r^2 - 2r + 0.5)/r^2); plot!(aspect_ratio=:equal)

function multistepcoefficients(m,n)
  s = length(m) + length(n) - 1
  A = (m.+1).^(0:s-1) .|> Rational
  B = (0:s-1).*((n.+1).^[0;0:s-2])
  c = -[A[:,2:end] B]\ones(Int64,s)
  [1;c[1:length(m)-1]], c[length(m):end] 
end

function plotstability(a,b)
  λk(r) = (a ⋅ r.^-(0:length(a)-1)) ./ (b ⋅ r.^-(0:length(b)-1))
  r = exp.(im*LinRange(0,2π,200))
  plot(λk.(r),label="",aspect_ratio=:equal)
end

m = [0 1]; n = [0 1 2]
a, b = zeros(maximum(m)+1), zeros(maximum(n)+1) 
a[m.+1], b[n.+1] = multistepcoefficients(m,n)
plotstability(a,b)

using DifferentialEquations, Plots
pendulum(u,p,t) = [u[2]; -sin(u[1])]
u₀ = [8π/9,0]; tspan = [0,8π]
problem = ODEProblem(pendulum, u₀, tspan)
method = Trapezoid()
solution = solve(problem,method)
plot(solution, xaxis="t", label=["θ" "ω"])

using DifferentialEquations, ParameterizedFunctions, DiffEqDevTools, Sundials
f = @ode_def LotkaVolterra begin
  dx = α*x - β*x*y
  dy = -γ*y + δ*x*y
end α β γ δ

params = [1.5,1.0,3.0,1.0]
problem = ODEProblem(f,[1.0;1.0],(0.0,10.0),params)

solution = solve(problem,Vern7(),abstol=1e-15,reltol=1e-15)
testsol = TestSolution(solution)
plot(solution)

solvers = [ 
  Dict(:alg=>VCABM())
  Dict(:alg=>DP5())
  Dict(:alg=>BS3())
  Dict(:alg=>CVODE_Adams())
  Dict(:alg=>Tsit5())
  Dict(:alg=>DP8())
  Dict(:alg=>Vern8())
  Dict(:alg=>Vern9())
];

abstols = 10 .^[-(LinRange(6,13,40))...]
reltols = 1e3*abstols
wps = WorkPrecisionSet(problem,abstols,reltols,solvers;
  appxsol=testsol,save_everystep=true,numruns=100,maxiters=10^5)

using DataFrames, LinRegOutliers
function logfit(w, c = 1)
  df = DataFrame(x=log.(w.errors), y=log.(w.times))
  reg = createRegressionSetting(@formula(y ~ x), df)
  s = lts(reg; crit = c)
  β = s["betas"]
  delete!(df, s["outliers"])
  x = exp.([minimum(df.x) maximum(df.x)])  
  y = [exp(β[1])*(x[1])^β[2] exp(β[1])*(x[2])^β[2]]  
  (x,y,exp.(df.x),exp.(df.y))
end

plot(xaxis=:log,yaxis=:log,xflip = true)
for i∈1:length(wps)
  x, y, xₚ, yₚ = logfit(wps[i],5)
  scatter!(xₚ,yₚ,alpha=0.25,labels=:none,markersize=5,color=i)
  plot!(x',y',linewidth=3,labels=:none,color=i)
  annotate!(x[1], y[1], text.(wps.names[i],9,:left))  
end
plot!()

using DifferentialEquations
θ₀ = π/3; ℓ = 1; tspan = (0,30.0)
U₀ = [ℓ*sin(θ₀), -ℓ*cos(θ₀), 0, 0]
function pendulum(dU,U,p,t)
  x,y,u,v = U; ℓ,g,m = p
  τ = m*(u^2 + v^2 + g*y)/ℓ^2
  dU .= (u, v, -τ*x/m, -τ*y/m + g)
end
problem = ODEProblem(pendulum, U₀, tspan, (ℓ,1,1))
solution = solve(problem, Tsit5())
plot(solution, idxs=(1,2), yflip=true, aspect_ratio=:equal)

function residual(r,w,p,t)
  x,y,u,v = w;  ℓ,g,m = p
  r .= (x^2+y^2-ℓ^2, x*u+y*v, 0, 0)
end
cb = ManifoldProjection(residual)
solution = solve(problem, Tsit5(), callback=cb)
plot(solution, idxs=(1,2), yflip=true, aspect_ratio=:equal)

using DifferentialEquations
θ₀ = π/3; ℓ = 1; tspan = (0,30.0)
u₀ = [ℓ*sin(θ₀), -ℓ*cos(θ₀), 0, 0, -ℓ*cos(θ₀)]
function pendulum(dU,U,p,t)
  x,y,u,v,τ = U; ℓ,g,m = p
  dU .= (u, v, -τ*x/m, -τ*y/m + g, -ℓ^2*τ + m*g*y + m*(u^2 + v^2))
end
M = Diagonal([1,1,1,1,0])
f = ODEFunction(pendulum, mass_matrix=M)
problem = ODEProblem(f,u₀,tspan,(ℓ,1,1))
solution = solve(problem, Rodas5(), reltol=1e-8, abstol=1e-8)
plot(solution, idxs=(1,2), yflip=true, aspect_ratio=:equal)

using ModelingToolkit, DifferentialEquations
θ₀ = π/3; ℓ = 1; tspan = (0,30.0)
u₀ = [ℓ*sin(θ₀), -ℓ*cos(θ₀), 0, 0, 0]
M = Diagonal([1,1,1,1,0])
function pendulum(dU,U,p,t)
  x,y,u,v,τ = U; ℓ,g,m = p
  dU .= (u, v, -τ*x/m, -τ*y/m + g, x^2 + y^2 - ℓ^2)
end
f = ODEFunction(pendulum, mass_matrix=M)
problem = ODEProblem(f, u₀, tspan, [ℓ,1,1])
sys = modelingtoolkitize(problem)
pendulum_sys = structural_simplify(dae_index_lowering(sys))
problem = ODAEProblem(pendulum_sys, [], tspan)
solution = solve(problem, Tsit5(), abstol=1e-8, reltol=1e-8)
plot(solution, idxs=(1,3), yflip=true, aspect_ratio=:equal)

using FastGaussQuadrature
function differentiation_matrix(n,Δt=1) 
  nodes, _ = gausslobatto(n+1)
  t = (nodes[2:end].+1)/2
  A = t.^(0:n-1)'.*(1:n)'
  B = t.^(1:n)'
  (A/B)/Δt, t
end

function pendulumGLL(U,(U₀,D,n))
  x,y,u,v,τ = U[1:n],U[n+1:2n],U[2n+1:3n],U[3n+1:4n],U[4n+1:5n]
  x₀,y₀,u₀,v₀,_ = U₀
  ℓ,g,m = (1,1,1)
  [D(x,x₀) .- u;
  D(y,y₀) .- v;
  D(u,u₀) .+ τ.*x/m;
  D(v,v₀) .+ τ.*y/m .- g;
  x.^2 + y.^2 .- ℓ^2]
end

using NonlinearSolve
θ₀ = π/3; ℓ = 1; tspan = (0,30.0)
u₀ = [ℓ*sin(θ₀), -ℓ*cos(θ₀), 0, 0, 0]
n = 5; N = 100; Δt = 30/N
M,_ = differentiation_matrix(n,Δt) 
D = (u,u₀) -> M*(u .- u₀)
u = u₀
for i = 2:N
  problem = NonlinearProblem(pendulumGLL,[ones(n)*u₀'...],(u₀,D,n))
  solution = solve(problem,NewtonRaphson(),abstol=1e-12)
  u = [u reshape(solution.u,5,n)']  
  u₀ = u[:,end]
end
plot(u[1,:], u[2,:], yflip=true, aspect_ratio=:equal)

Δx = .01; Δt = .01; L = 2; ν = Δt/Δx^2; uₗ = 0; uᵣ = 0;
x = -L:Δx:L; m = length(x) 
u = (abs.(x).<1)
u[1] += 2ν*uₗ; u[m] += 2ν*uᵣ 
D = Tridiagonal(ones(m-1), -2ones(m), ones(m-1))
D[1,2] = 0; D[m,m-1] = 0
A = I - ν*D 
for i in 1:20
  u = A\u
end
plot(x,u)

Δx = .01; Δt = .01; L = 2; ν = Δt/Δx^2
x = -L:Δx:L; m = length(x) 
u = float.(abs.(x).<1)
D = Tridiagonal(ones(m-1), -2ones(m), ones(m-1))
D[1,2] = 2; D[m,m-1] = 2
A = 2I + ν*D 
B = 2I - ν*D 
anim = @animate for i in 1:40
  plot(x,u, legend=:none, ylims=(0,1))
  u .= B\(A*u)
end
gif(anim, "heat_cn.gif", fps = 5)

using Sundials
m = 400; L = 2; x = LinRange(-L,L,m); Δx = x[2]-x[1]
α = (u -> u.^2)
Du(u,Δx,t) = [0;diff(α((u[1:end-1]+u[2:end])/2).*diff(u))/Δx^2;0]
u₀ = (abs.(x).<1)
problem = ODEProblem(Du,u₀,(0,2),Δx)
method = CVODE_BDF(linear_solver=:Band, jac_lower=1, jac_upper=1)
solution = solve(problem, method);

anim = @animate for t in LinRange(0,2,200)
  plot(x,solution(t),legend=:none,fill=(0,0.4,:red),ylims=(0,1))
end
gif(anim, "porous.gif", fps = 15)

using Interact
@manipulate for t∈slider(0:0.01:2; value=0, label="time")
  plot(x,solution(t),  fill = (0, 0.4, :red))
  plot!(ylims=(0,1),legend=:none)
end

using FFTW
m = 256; ℓ = 4
ξ² = (fftshift(-(m-1)÷2:m÷2)*(2π/ℓ)).^2
u = (t,u₀) -> real(ifft(exp.(-ξ²*t).*fft(u₀)))

x = LinRange(-ℓ/2,ℓ/2,m+1)[2:end]
u₀ = (abs.(x).<1)
plot(x,u(0.1,u₀))

F = plan_fft(u₀)
u = (t,u₀) -> real(F\(exp.(-ξ²*t).*(F*u₀)))

#= In moving from Julia 1.7 to 1.9 this code no longer works.
using FFTW, DifferentialEquations, SciMLOperators
ℓ = 100; T = 200.0; n = 512
x = LinRange(-ℓ/2,ℓ/2,n+1)[2:end]
u₀ = exp.(-x.^2)
iξ = im*[0:(n÷2);(-n÷2+1):-1]*2π/ℓ
F = plan_fft(u₀)
L = DiagonalOperator(-iξ.^2-iξ.^4)
N = (u,_,_) -> -0.5iξ.*(F*((F\u).^2))
problem = SplitODEProblem(L,N,F*u₀,(0,T))
method = KenCarp3(linsolve=KrylovJL(),autodiff=false)
solution = solve(problem,method,reltol=1e-12)
u = t -> real(F\solution(t))
=#

#=
using ImageShow,ColorSchemes
s = hcat(u.(LinRange(0,T,500))...)
get(colorschemes[:binary],abs.(s),:extrema)
=#

Δᵒ(Q,step=1) = Q - circshift(Q,(step,0))
flux(Q,c) = c.*Δᵒ(Q) - 0.5c.*(1 .- c).*(Δᵒ(Q,1)+Δᵒ(Q,-1))
H = (u,v,iξ₁,iξ₂) -> F*((F\u).*(F\(iξ₁.*u)) + (F\v).*(F\(iξ₂.*u)))

using FFTW
ℓ = 2; n = 128; ϵ = 0.001; Δt = 0.001; Δx = ℓ/n
x = LinRange(Δx,ℓ,n)'; y = x'
Q = (@. 0.5(1 + tanh(10 - 20abs(1 - 2y/ℓ)))).*ones(1,n)
uⁿ = Q .* (1 .+ 0.5sin.(ℓ*π*x))
vⁿ = zeros(n,n) 
F = plan_fft(uⁿ)
uⁿ, vⁿ = F*uⁿ, F*vⁿ 
uᵒ, vᵒ = uⁿ, vⁿ
iξ₁ = im*fftfreq(n,n)'*(2π/ℓ)
iξ₂ = transpose(iξ₁)
ξ² = iξ₁.^2 .+ iξ₂.^2
H₁ⁿ, H₂ⁿ = H(uⁿ,vⁿ,iξ₁,iξ₂),  H(vⁿ,uⁿ,iξ₂,iξ₁)
M₊ = (1/Δt .+ (ϵ/2)*ξ²)
M₋ = (1/Δt .- (ϵ/2)*ξ²);

for i = 1:1200
  Q -= flux(Q,(Δt/Δx).*real(F\vⁿ)) + flux(Q',(Δt/Δx).*real(F\uⁿ)')'
  H₁ⁿ⁻¹, H₂ⁿ⁻¹ = H₁ⁿ, H₂ⁿ 
  H₁ⁿ, H₂ⁿ = H(uⁿ,vⁿ,iξ₁,iξ₂),  H(vⁿ,uⁿ,iξ₂,iξ₁)
  uᵒ = uⁿ - uᵒ + (-1.5H₁ⁿ + 0.5H₁ⁿ⁻¹ + M₊.*uⁿ)./M₋
  vᵒ = vⁿ - vᵒ + (-1.5H₂ⁿ + 0.5H₂ⁿ⁻¹ + M₊.*vⁿ)./M₋
  ϕ = (iξ₁.*uᵒ + iξ₂.*vᵒ)./(ξ² .+ (ξ².≈ 0)) 
  uⁿ, vⁿ = uᵒ - iξ₁.*ϕ, vᵒ - iξ₂.*ϕ
end

using ColorSchemes
imresize(get(ColorSchemes.acton, clamp01.(Q)),ratio=3)

using LinearAlgebra, Plots
N = 10000
n = [sum([!(det(rand(Bool,d,d))≈0) 
 for i in 1:N]) for d in 1:20]
plot(n/N,marker=:o)
p = [1,6,174,22560,12514320,28836612000,270345669985440,10160459763342013440]
r = @. 2^(log2(p)-(1:8)^2)
plot!(r,marker=:square,opacity=0.5)

A = {{0,1,0,1},{1,0,2,1},{0,2,0,0},{1,1,0,0}};
W = Inverse[IdentityMatrix[4] - z*A];
Simplify[W[[1,3]]]
Series[%,{z,0,10}]

A = [0 1 0 1;1 0 2 1;0 2 0 0;1 1 0 0]
walks(A,i,j,N) = [(A^n)[i,j] for n in 1:N]
walks(A,1,3,9)

function detx(A)
  L,U,P = lu(A)
  s = 1
  for i in 1:length(P)
    m = findfirst(P.==i)
    i!=m && ( s *= -1; P[[m,i]]=P[[i,m]] )
  end
  s * prod(diag(U))
end

A = rand(20,20)
detx(A) ≈ det(A)

function rcuthillmckee(A)
  r = sortperm(vec(sum(A.!=0,dims=2)))
  p = Int64[]
  while ~isempty(r)
    q = [popfirst!(r)]
    while ~isempty(q)
      q₁ = popfirst!(q)
      append!(p,q₁)  
      k = findall(!iszero, A[q₁,r])
      append!(q,splice!(r,k))
    end
  end
  reverse(p)
end

using SparseArrays
A = sparse(Symmetric(sprand(1000,1000,0.003)))
p = rcuthillmckee(A)
plot(spy(A), spy(A[p,p]), colorbar = false)

using DelimitedFiles, SparseArrays, Graphs, GraphPlot
ij = readdlm(download(bucket*"dolphins.csv"),',',Int)
A = sparse(ij[:,1],ij[:,2],one.(ij[:,1]))
gplot(Graph(A),layout=circular_layout)

p = rcuthillmckee(A)
gplot(Graph(A[p,p]),layout=circular_layout)

using CSV, DataFrames
diet = DataFrame(CSV.File(download(bucket*"diet.csv")))

A = Array(diet[2:end,4:end])'
c = ones.(size(A,2))
b = Array(diet[1,4:end])
food = diet[2:end,1]
solution = simplex(b,A',c)
print("foods: ", food[solution.y .!= 0], "\n")
print("daily cost: ", solution.z)

using JuMP, GLPK
model = Model(GLPK.Optimizer)
@variable(model, x[1:size(A,2)] ≥ 0)
@objective(model, Min, c' * x)
@constraint(model, A * x .≥ b)
optimize!(model)
print("foods: ", food[JuMP.value.(x) .!= 0], "\n")
print("daily cost: ", objective_value(model))

function findpath(A,a,b)
  r = collect(1:size(A,2))
  q = [a]; p = [-1]; i = 0
  splice!(r,a)
  while i<length(q)
    k = findall(!iszero, A[q[i+=1],r])
    any(r[k].==b) && return(append!(q[backtrack(p,i)],b))
    append!(q,splice!(r,k))
    append!(p,fill(i,length(k)))         
  end     
  display("No path.")  
end

function backtrack(p,i)  
  s = [];  while (i!=-1); append!(s,i); i = p[i]; end
  return(reverse(s))
end

using DelimitedFiles, SparseArrays
get_names(filename) = readdlm(download(bucket*filename*".txt"),'\n')
function get_adjacency_matrix(filename)
  i = readdlm(download(bucket*filename*".csv"),',',Int,'\n')
  sparse(i[:,1], i[:,2], 1)
end

actors = get_names("actors")
movies = get_names("movies")
B = get_adjacency_matrix("actor-movie");

(m,n) = size(B)
A = [spzeros(n,n) B' ; B spzeros(m,m)];
actormovie = [ actors ; movies ];

index = x -> findfirst(actors.==x)[1]
actormovie[findpath(A, index("Emma Watson"), index("Kevin Bacon"))]

actormovie[findpath(A, index("Bruce Lee"), index("Tommy Wiseau"))]

d = 0.5; A = sprand(800,600,d)
Base.summarysize(A)/Base.summarysize(Matrix(A))

X = load(download(bucket*"laura.png")) .|> Gray
m,n = size(X)
blur = x -> exp(-(x/20).^2/2)
A = [blur(i-j) for i=1:m, j=1:m]; A ./= sum(A,dims=2)
B = [blur(i-j) for i=1:n, j=1:n]; B ./= sum(B,dims=1)
N = 0.01*rand(m,n)
Y = A*X*B + N;

α = 0.05
X₁ = A\Y/B
X₂ =  (A'*A+α^2*I)\A'*Y*B'/(B*B'+α^2*I)
X₃ = pinv(A,α)*Y*pinv(B,α)
Gray.([X Y X₁ X₂ X₃])

vandermonde(t,n) = vec(t).^(0:n-1)'
build_poly(c,X) = vandermonde(X,length(c))*c

function solve_filip(x,y,n)
  V = vandermonde(x, n)
  c = Array{Float64}(undef, 3, n)
  c[1,:] = (V'*V)\(V'*y)
  Q,R = qr(V)
  c[2,:] = R\(Matrix(Q)'*y)
  c[3,:] = pinv(V,1e-9)*y
  r = [norm(V*c[i,:]-y) for i in 1:3]
  return(c,r)
end

using DelimitedFiles
data = readdlm(download(bucket*"filip.csv"),',',Float64)
coef = readdlm(download(bucket*"filip-coeffs.csv"),',')
(x,y) = (data[:,2],data[:,1])
β,r = solve_filip(x, y, 11)
X = LinRange(minimum(x),maximum(x),200)
Y = [build_poly(β[i,:],X) for i in 1:3]
plot(X,Y); scatter!(x,y,opacity=0.5)
[coef β']

using Statistics
zscore(X,x=X) = (X .- mean(x))/std(x)
c,r = solve_filip(zscore(x), zscore(y), 11)
Y = [build_poly(c[i,:],zscore(X,x))*std(y).+mean(y) for i in 1:3]
plot(X,Y); scatter!(x,y,opacity=0.5)

using CSV, DataFrames, Dates
data = download(bucket*"dailytemps.csv")
df = DataFrame(CSV.File(data))
day = Dates.value.(df[:,:date] .- df[1,:date])/365
u = df[:,:temperature]
tempsmodel(t) = [sin.(2π*t) cos.(2π*t) one.(t)]
c = tempsmodel(day)\u
scatter(df[:,:date],u,alpha=0.3)
plot!(df[:,:date],tempsmodel(day)*c,width=3)

using MLDatasets, Arpack, Images
image_train, label_train = MNIST(split=:train)[:]
image_train = reshape(permutedims(image_train,[3,2,1]),60000,:)
V = [( D = image_train[label_train.==i,:];
  svds(D,nsv=12)[1].V ) for i ∈ 0:9];

pix = -abs.(reshape(V[3+1],28,:))
rescale = scaleminmax(minimum(pix), maximum(pix))
pix .|> rescale .|> Gray

image_test, label_test = MNIST(split=:test)[:]
image_test = reshape(permutedims(image_test,[3,2,1]),10000,:)
r = [( q = image_test*V[i]*V[i]' .- image_test;
  sum(q.^2,dims=2) ) for i∈1:10]
r = hcat(r...) 
prediction = [argmin(r[i,:]) for i∈1:10000] .- 1
[sum(prediction[label_test.== i].== j) for j∈0:9, i∈0:9]

actors = get_names("actors")
genres = get_names("genres")
A = get_adjacency_matrix("movie-genre")
A = (A*Diagonal(1 ./ sum(A,dims=1)[:]))
B = get_adjacency_matrix("actor-movie");

using Arpack
X,_ = svds(A*B,nsv=10)
U = X.U; Σ = Diagonal(X.S); Vᵀ = X.Vt
Q = Vᵀ./sqrt.(sum(Vᵀ.^2,dims=1));

index = x -> findfirst(actors.==x)[1]
q = Q[:,index("Steve Martin")][:]
z = Q'*q
r = sortperm(-z)
actors[r[1:10]]

p = U*Σ*q
r = sortperm(-p)
[genres[r] p[r]/sum(p)]

X = [3 3 12; 1 15 14; 10 2 13; 12 15 14; 0 11 12]
reference = X[1,:]; X .-= reference'
A = [2 2 -2].*X; b = (X.^2)*[1; 1; -1]
xₒₗₛ = A\b + reference
xₜₗₛ = tls(A,b) + reference
xₒₗₛ, xₜₗₛ

using CSV, DataFrames
df = DataFrame(CSV.File(download(bucket*"shotspotter.csv")))
X = Array(df[1:end-1,:]); xₒ = Array(df[end,:]);

c² = 328.9^2
reference = X[1,:]; X .-= reference'
A = [2 2 2 -2c²].*X; b = (X.^2)*[1; 1; 1; -c²]
xₒₗₛ = A\b + reference
error = xₒₗₛ - xₒ

n = 8
E = collect(Iterators.flatten([eigvals(randn(n,n)) for i∈1:2000]))
scatter(E,mc=:black,ma=0.05,legend=nothing,aspect_ratio=:equal)

function rayleigh(A)
  x = randn(size(A,1),1)
  while true
    x = x/norm(x) 
    ρ = (x'*A*x)[1]
    M = A - ρ*I
    abs(cond(M, 1)) < 1e12 ? x = M\x : return(ρ,x)
  end
end
A = [2 3 6 4; 3 0 3 1; 6 3 8 8; 4 1 8 2]
rayleigh(A)

function implicitqr(A)
  n = size(A,1)
  tolerance = 1E-12
  H = Matrix(hessenberg(A).H)
  while true
    if abs(H[n,n-1])<tolerance
      if (n-=1)<2;  return(diag(H)); end
    end
    Q = givens([H[1,1]-H[n,n];H[2,1]],1,2)[1]
    H[1:2,1:n] = Q*H[1:2,1:n]
    H[1:n,1:2] = H[1:n,1:2]*Q'
    for i = 2:n-1
      Q = givens([H[i,i-1];H[i+1,i-1]],1,2)[1]
      H[i:i+1,1:n] = Q*H[i:i+1,1:n]
      H[1:n,i:i+1] = H[1:n,i:i+1]*Q'
    end
  end
end
n = 20; S = randn(n,n)
D = Diagonal(1:n); A = S*D*inv(S)
implicitqr(A)

actors = get_names("actors")
B = get_adjacency_matrix("actor-movie")
M = (B'*B .!= 0) - I
v = ones(size(M,1))
for k = 1:8
  v = M*v; v /= norm(v);
end
r = sortperm(-v); actors[r][1:10]

findfirst(actors[r] .== "Kevin Bacon")

actors[ -sum(M,dims=1)[:] |> sortperm ][1:10]

function randomizedsvd(A,k)
  Z = rand(size(A,2),k);
  Q = Matrix(qr(A*Z).Q)
  for i in  1:3
    Q = Matrix(qr(A'*Q).Q)
    Q = Matrix(qr(A*Q).Q)
  end
  W = svd(Q'*A)
  return((Q*W.U,W.S,W.V))  
end

using Images
img = load(download(bucket*"red-fox.jpg"))
A = Float64.(Gray.(img))
U,S,V = randomizedsvd(A,10)
Gray.([A U*Diagonal(S)*V'])

using Arpack
@time randomizedsvd(A,10)
@time svds(A, nsv=10)
@time svd(A);

ϵ =  1 .- sqrt.(cumsum(S.^2))/norm(A)
plot(ϵ,marker=:circle)

using Arpack
m = 5000
A = [1/(1 + (i+j-1)*(i+j)/2 - j) for i∈1:m, j∈1:m]
svds(A,nsv=1)[1].S[1]

using SparseArrays, LinearAlgebra
⊗(x,y) = kron(x,y); ϕ = x -> x-x^2
n = 50 ; x = (1:n)/(n+1); Δx = 1/(n+1)
J = sparse(I, n, n)
D = spdiagm(-1 => ones(n-1), 0 => -2ones(n), 1 => ones(n-1) )
A = ( D⊗J⊗J + J⊗D⊗J + J⊗J⊗D ) / Δx^2
f = [ϕ(x)*ϕ(y) + ϕ(x)*ϕ(z) + ϕ(y)*ϕ(z) for x∈x, y∈x, z∈x][:]
uₑ = [ϕ(x)*ϕ(y)*ϕ(z)/2 for x∈x, y∈x, z∈x][:];

function stationary(A,b,ω=0,n=400)
  ϵ = []; u = zero.(b)
  P = (1-ω)*sparse(Diagonal(A)) + ω*sparse(LowerTriangular(A))
  for i=1:n
    u += P\(b-A*u)
    append!(ϵ,norm(u - uₑ,1))
  end
  return(ϵ)
end

function conjugategradient(A,b,n=400)
  ϵ = []; u = zero.(b)  
  p = r = b-A*u
  for i=1:n
    Ap = A*p
    α = (r'*p)/(Ap'*p)
    u += α.*p; r -= α.*Ap
    β = (r'*Ap)/(Ap'*p)
    p = r - β.*p
    append!(ϵ,norm(u - uₑ,1))
  end
  return(ϵ)
end

ϵ = zeros(400,4)
@time ϵ[:,1] = stationary(A,-f,0)
@time ϵ[:,2] = stationary(A,-f,1)
@time ϵ[:,3] = stationary(A,-f,1.9)
@time ϵ[:,4] = conjugategradient(A,-f)
method = ["Jacobi" "Gauss-Seidel" "SOR" "Conjugate Gradient"]
plot(ϵ,yaxis=:log,labels=method,bglegend=RGBA(1,1,1,0.7))

k = 1:120;  ([one.(k) k]\log10.(ϵ[k,:]))[2,:]

using Primes, LinearAlgebra, SparseArrays, IterativeSolvers
n = 20000
d = 2 .^ (0:14); d = [-d;d]
P = Diagonal(primes(224737))
B = [ones(n - abs(d)) for d∈d]
A = sparse(P) + spdiagm(n ,n, (d .=> B)...)
b = zeros(n); b[1] = 1
cg(A, b; Pl=P)[1]

function fftx3(c)
  n = length(c)
  ω = exp(-2im*π/n)
  if mod(n,3) == 0 
    k = collect(0:n/3-1)
    u = [transpose(fftx3(c[1:3:n-2]));
         transpose((ω.^k).*fftx3(c[2:3:n-1]));
         transpose((ω.^2k).*fftx3(c[3:3:n]))]  
    F = exp(-2im*π/3).^([0;1;2]*[0;1;2]')
    return(reshape(transpose(F*u),:,1))
  else 
    F =  ω.^(collect(0:n-1)*collect(0:n-1)')
    return(F*c)
  end 
end

using FFTW
v = rand(24,1)
[fft(v) fftx3(v)]

using FFTW
function multiply(p_,q_)
  p = [parse.(Int,split(reverse(p_),""));zeros(length(q_),1)]
  q = [parse.(Int,split(reverse(q_),""));zeros(length(p_),1)]
  pq = Int.(round.(real.(ifft(fft(p).*fft(q)))))
  carry = pq .÷ 10
  while any(carry.>0)
    pq -= carry*10 - [0;carry[1:end-1]]
    carry = pq .÷ 10
  end
  n = findlast(x->x>0, pq)
  return(reverse(join(pq[1:n[1]])))
end

p = "32769132993266709549961988190834461413177642967992942539798288533"
q = "3490529510847650949147849619903898133417764638493387843990820577"
multiply(p,q), parse(BigInt, p)*parse(BigInt, q)

using Random
p = randstring('0':'9', 100000)
q = randstring('0':'9', 100000)
@time multiply(p,q)
@time parse(BigInt, p)*parse(BigInt, q);

function dctII(f)
  n = size(f,1)
  ω = exp.(-0.5im*π*(0:n-1)/n)
  return real(ω.*fft(f[[1:2:n; n-mod(n,2):-2:2],:],1))
end

function idctII(f)
  n = size(f,1)
  ω = exp.(-0.5im*π*(0:n-1)/n)
  f[1,:] = f[1,:]/2
  f[[1:2:n; n-mod(n,2):-2:2],:] = 2*real(ifft(f./ω,1))
  return f
end

dct2(f) = dctII(dctII(f')')
idct2(f) = idctII(idctII(f')')

function dctcompress2(A,d)
  B = dct(Float64.(A))
  m,n = size(A)  
  m₀,n₀ = sqrt(d).*(m,n) .|> floor .|> Int
  B₀ = B[1:m₀,1:n₀]
  A₀ = idct([B₀ zeros(m₀,n-n₀); zeros(m-m₀,n)])
  A₀ = A₀ |> clamp01! .|> Gray
  sizeof(B₀)/sizeof(B), sqrt(1-(norm(B₀)/norm(B))^2), A₀
end

A = Gray.(load(download(bucket*"laura_square.png")))
[A dctcompress2(A,0.05)[3]]

aitken₁(x₁,x₂,x₃) = x₃ - (x₃-x₂)^2/(x₃ - 2x₂ + x₁)
aitken₂(x₁,x₂,x₃) = (x₁*x₃ - x₂^2)/(x₃ - 2x₂ + x₁)

n = 60000
p = cumsum([(-1)^i*4/(2i+1) for i=0:n])
p₁ = aitken₁.(p[1:n-2],p[2:n-1],p[3:n])
p₂ = aitken₂.(p[1:n-2],p[2:n-1],p[3:n])
plot(abs.(π.-p)); plot!(abs.(π.-p₂)); plot!(abs.(π.-p₁))
plot!(xaxis=:log,yaxis=:log)

using DifferentialEquations
function homotopy(f,df,x)
  dxdt(x,p,t) = -df(x)\p
  sol = solve(ODEProblem(dxdt,x,(0.0,1.0),f(x)))
  sol.u[end]
end

function newton(f,df,x)
  for i in 1:100
    Δx = -df(x)\f(x)
    norm(Δx) > 1e-8 ? x += Δx : return(x)
  end
end

f = z -> ((x,y)=tuple(z...);
  [(x^2+y^2)^2-2(x^2-y^2); (x^2+y^2-1)^3-x^2*y^3])
df = z -> ((x,y)=tuple(z...);
  [4x*(x^2+y^2-1)  4y*(x^2+y^2+1);
   6x*(x^2+y^2-1)^2-2x*y^3  6y*(x^2+y^2-1)^2-3x^2*y^2])

display(homotopy(f,df,[1,1]))
display(newton(f,df,[1,1]))

function ⊕(P,Q)
  a = 0
  r = BigInt(2)^256 - 2^32 - 977
  if P[1]==Q[1]
    d = invmod(2*P[2], r)
    λ = mod((3*powermod(P[1],2,r)+a)*d, r)
  else
    d = invmod(Q[1]-P[1], r)
    λ = mod((Q[2]-P[2])*d, r)
  end
  x = mod(powermod(λ, 2, r) - P[1] - Q[1], r)
  y = mod(-λ*(x-P[1])-P[2], r)
  [x;y]
end

Base.isodd(m::BigInt) = ((m&1)==1)
function dbl_add(m,P)
  if m > 1
    Q = dbl_add(m>>1,P)
    return isodd(m) ? (Q⊕Q)⊕P : Q⊕Q
  else
    return P
  end
end

P₁ = big"0x79BE667EF9DCBBAC55A06295CE87
        0B07029BFCDB2DCE28D959F2815B16F81798"
P₂ = big"0x483ADA7726A3C4655DA4FBFC0E11
        08A8FD17B448A68554199C47D08FFB10D4B8"
P = [P₁; P₂]
m, n = 1829442, 3727472
mP = dbl_add(m,P)
nmP = dbl_add(n,mP)
nP = dbl_add(n,P)
print("Alice's private key:\n$m\n\n")
print("Bob's private key:\n$n\n\n")
print("Alice's public key:\n$mP\n\n")
print("Bob's public key:\n$nP\n\n")
print("Alice and Bob's shared secret key:\n$nmP")

function spline_periodic(x,y)
  h = diff(x) 
  d = 6*diff(diff([y[end-1];y])./[h[end];h])
  α = h[1:end-1]
  β = 2*(h+circshift(h,1))
  C = Matrix(SymTridiagonal(β,α))
  C[1,end]=h[end]; C[end,1]=h[end] 
  m = C\d
  return([m;m[1]])
end

function evaluate_spline(xᵢ,yᵢ,m,n)
  h = diff(xᵢ)
  B = yᵢ[1:end-1] .- m[1:end-1].*h.^2/6
  A = diff(yᵢ)./h - h./6 .*diff(m)
  x = range(minimum(xᵢ),maximum(xᵢ),length=n+1)
  i = sum(xᵢ' .≤ x,dims=2)
  i[end] = length(xᵢ)-1 
  y = @. (m[i]*(xᵢ[i+1]-x)^3 + m[i+1]*(x-xᵢ[i])^3)/(6h[i]) +
     A[i]*(x-xᵢ[i]) + B[i]
  return(x,y)
end

n = 5; nx = 30
x = rand(n); y = rand(n)
x = [x;x[1]]; y = [y;y[1]]
t = [0;cumsum(sqrt.(diff(x).^2 + diff(y).^2))]
X = evaluate_spline(t,x,spline_periodic(t,x),nx*n)
Y = evaluate_spline(t,y,spline_periodic(t,y),nx*n)
scatter(x,y); plot!(X[2],Y[2],legend=false)

n = 20; N = 200
x = LinRange(-1,1,n)
y = float(x .> 0);

ϕ₁(x,a) = abs.(x.-a).^3
ϕ₂(x,a) = exp.(-20(x.-a).^2)
ϕ₃(x,a) = x.^a

X = LinRange(-1,1,N)
interp(ϕ,a)  = ϕ(X,a')*(ϕ(x,a')\y)
Y₁ = interp(ϕ₁,x)
Y₂ = interp(ϕ₂,x)
Y₃ = interp(ϕ₃,(0:n-1))
scatter(x,y,seriestype = :scatter, marker = :o, legend = :none)
plot!(X,[Y₁,Y₂,Y₃],ylim=[-0.5,1.5])

function collocation_solve(L,f,bc,x)
  h = x[2]-x[1]
  S = L(x)*([1 -1/2 1/6; -2 0 2/3; 1 1/2 1/6]./[h^2 h 1])'
  d = [bc[1]; f(x); bc[2]]
  A = Matrix(Tridiagonal([S[:,1];0], [0;S[:,2];0], [0;S[:,3]]))
  A[1,1:3], A[end,end-2:end] = [1 4 1]/6, [1 4 1]/6 
  return(A\d)
end

function collocation_build(c,x,N)
  X = LinRange(x[1],x[end],N)
  h = x[2] - x[1]
  i = Int32.(X .÷ h .+ 1); i[N] = i[N-1]
  C = [c[i] c[i.+1] c[i.+2] c[i.+3]]'
  B = (x->[(1-x).^3;4-3*(2-x).*x.^2;4-3*(1+x).*(1-x).^2;x.^3]/6)
  Y = sum(C.*hcat(B.((X.-x[i])/h)...),dims=1)
  return(Array(X),reshape(Y,(:,1)))
end

using Roots, SpecialFunctions
n = 20; N = 141
L = (x -> [x one.(x) x])
f = (x -> zero.(x) )
b = fzero(besselj0, 11)
x = range(0,b,length=n)
c = collocation_solve(L,f,[1,0],x)
X,Y = collocation_build(c,x,70)
plot(X,[Y besselj0.(X)])

N = 10*2 .^(1:7); ϵ = []
for n in N
  x = LinRange(0,b,n) 
  c = collocation_solve(L,f,[1,0],x)
  X,Y = collocation_build(c,x,n)
  append!(ϵ, norm(Y-besselj0.(X)) / n)
end
slope = ([log.(N) one.(N)]\log.(ϵ))[1]
s = "slope = "*string.(round.(slope,digits=2))
plot(N, ϵ, xaxis=:log, yaxis=:log, marker=:o, label=s)

using FFTW
n = 256; ℓ = 2
x = (0:n-1)/n*ℓ .- ℓ/2
ξ = [0:(n/2-1); (-n/2):-1]*(2π/ℓ)
f₁ = exp.(-16*x.^2)
f₂ = sin.(π*x)
f₃ = x.*(1 .- abs.(x))
deriv(f,p) = real(ifft((im*ξ).^p.*fft(f)))

using Interact
func = Dict("Gaussian"=>f₁,"sine"=>f₂,"spline"=>f₃)
@manipulate for f in togglebuttons(func; label="function"), 
  p in slider(0:0.01:2; value=0, label="derivative")
  plot(x,deriv(f,p),legend=:none)
end

using MLDatasets, Flux
model = Chain(
  Conv((5,5), 1=>6, tanh, pad=SamePad()),
  MeanPool((2,2)),
  Conv((5,5), 6=>16, tanh),
  MeanPool((2,2)),
  Conv((5,5), 16=>120, tanh),
  Flux.flatten,
  Dense(120, 84, tanh),
  Dense(84, 10))

image_train, label_train = MLDatasets.MNIST(split=:train)[:]
image_test, label_test = MLDatasets.MNIST(split=:test)[:]
image_train = Flux.unsqueeze(image_train, 3) 
image_test = Flux.unsqueeze(image_test, 3) 
label_train = Flux.onehotbatch(label_train, 0:9)
label_test = Flux.onehotbatch(label_test, 0:9);

loss(x,y) = Flux.Losses.logitcrossentropy(model(x),y)
parameters = Flux.params(model)
data = Flux.Data.DataLoader((image_train,label_train),batchsize=100)
optimizer = ADAM()

using ProgressMeter
@showprogress for epochs = 1:5
  Flux.train!(loss, parameters, data, optimizer)
end

accuracy(x,y) = sum(Flux.onecold(x) .== Flux.onecold(y))/size(y,2)
accuracy(model(image_test),label_test)

using LsqFit
ϵ = (x,p) -> @. √((x[:,1]-p[1])^2 + (x[:,2]-p[2])^2) - (x[:,3]-p[3])
x = [3 3 12; 1 15 14; 10 2 13; 12 15 14; 0 11 12]
curve_fit(ϵ,x,zeros(5),zeros(3)).param

using CSV, DataFrames
df = DataFrame(CSV.File(download(bucket*"shotspotter.csv")))
X = Array(df[1:end-1,:]); xₒ = Array(df[end,:]);

using LsqFit
c = 328.6
ϵ = (x,p)->sqrt.(sum([(x[:,i].-p[i]).^2 for i∈1:3]))-c*(x[:,4].-p[4])
xₙₗₛ = curve_fit(ϵ,X,zero(X[:,1]),X[1,:]).param
error = xₙₗₛ - xₒ

function plot_circle(x,y,r)
  t = LinRange(0,2π,100);
  plot!(x.+r*cos.(t), y.+r*sin.(t),color="black",opacity=0.5);
end

r = 328.6*(X[:,end] .- xₙₗₛ[end])
scatter(X[:,1], X[:,2], color = "black")
[plot_circle(X[i,1],X[i,2],r[i]) for i in 1:length(r)]
scatter!([xₙₗₛ[1]], [xₙₗₛ[2]], color = "red")
plot!(legend=:none,aspect_ratio=:equal)

d = [0,1,2,3]; n = length(d)
C = inv( d.^(0:n-1)' .// factorial.(0:n-1)' )
[C C*d.^n/factorial(n)]

function richardson(f,x,m)
  D = []
  for i in 1:m
    append!(D, ϕ(f,x,2^i))
    for j in i-1:-1:1 
      D[j] = (4^(i-j)*D[j+1] - D[j])/(4^(i-j) - 1)
    end
  end
  D[1]
end

ϕ = (f,x,n) -> (f(x+1/n) - f(x-1/n))/(2/n)
richardson(x->sin(x),0,4)

struct Dual
  value
  deriv
end

Dual(x) = Dual(x,1)
value(x) = x isa Dual ? x.value : x
deriv(x) = x isa Dual ? x.deriv : 0

Base.:+(u, v) = Dual(value(u)+value(v), deriv(u)+deriv(v))
Base.:-(u, v) = Dual(value(u)-value(v), deriv(u)-deriv(v))
Base.:*(u, v) = Dual(value(u)*value(v), 
   value(u)*deriv(v)+value(v)*deriv(u))
Base.:sin(u)  =  Dual(sin(value(u)), cos(value(u))*deriv(u))

Base.:sqrt(u) = Dual(sqrt(value(u)), deriv(u) / 2sqrt(value(u)))
Base.:/(u, v) = Dual(value(u)/value(v), 
  (value(u)*deriv(v)-value(v)*deriv(u))/(value(v)*value(v)))
Base.:cos(u)  = Dual(cos(value(u)), -sin(value(u))*deriv(u))

function get_zero(f,x)
  ϵ = 1e-12; δ = 1
  while abs(δ) > ϵ
    fx = f(Dual(x))
    δ = value(fx)/deriv(fx)
    x -= δ
  end
  return(x)
end

get_zero(x->4sin(x)+sqrt(x),4)

function get_extremum(f,x)
  ϵ = 1e-12; δ = 1
  while abs(δ)>ϵ
    fx = f(Dual(Dual(x)))
    δ = deriv(value(fx))/deriv(deriv(fx))
    x -= δ
  end
  return(x)
end

get_extremum(x->4sin(x)+sqrt(x),4)

function cauchyderivative(f, a, p, n = 20, ϵ = 0.1)
  ω = exp.(2π*im*(0:(n-1))/n)
  factorial(p)/(n*ϵ^p)*sum(@. f(a+ϵ*ω)/ω^p)
end

f = x -> exp(x)/(cos(x)^3 + sin(x)^3)
cauchyderivative(f, 0, 6)

function gauss_legendre(n)
  x = -cos.((4*(1:n).-1)*π/(4n+2))
  Δ = one.(x)
  dPₙ = 0
  while(maximum(abs.(Δ))>1e-16)
    Pₙ, Pₙ₋₁ = x, one.(x)
    for k ∈ 2:n
      Pₙ, Pₙ₋₁ = ((2k - 1)*x.*Pₙ-(k-1)*Pₙ₋₁)/k, Pₙ 
    end
    dPₙ = @. n*(x*Pₙ - Pₙ₋₁)/(x^2-1)
    Δ =  Pₙ ./ dPₙ 
    x -= Δ
  end
  return(x, @. 2/((1-x^2)*dPₙ^2))
end

f = (x-> 2sqrt(1-x^2))
x,w = gauss_legendre(10)
w ⋅ f.(x)

using FastGaussQuadrature
ξ,w = gausshermite(40)
u₀ = x ->  sin(x)
u = (t,x) -> w ⋅ u₀.(x.-2sqrt(t)*ξ)/sqrt(π)
x = LinRange(-12,12,100); plot(x,u.(1,x))

mc_π(n) = sum(sum(rand(2,n).^2,dims=1).<1)/n*4

m = 20; d = []; N = 2 .^ (1:20)
for n ∈ N
  append!(d,sum([abs(π - mc_π(n)) for i∈1:m])/m)
end
s = log.(N.^[0 1])\log.(d)
scatter(N,d,xaxis=:log, yaxis=:log) 
plot!(N,exp(s[1]).*N.^s[2])

using FastGaussQuadrature
function differentiation_matrix(n) 
  nodes, _ = gausslobatto(n+1)
  t = (nodes[2:end].+1)/2
  A = t.^(0:n-1)'.*(1:n)'
  B = t.^(1:n)'
  A/B,t
end

n = 20
M,t = differentiation_matrix(n) 
D = (u,u₀) -> M*(u .- u₀)

using NLsolve
u₀ = 1.0; α = 0.9
F = (u,u₀,α) -> D(u,u₀) .- α*u.^2
u = nlsolve(u->F(u,u₀,α),u₀*ones(n)).zero
plot([0;t], [u₀;u], marker=:circle, legend=false)

N = 30; Δt = 1.0/N; n = 8
M,t = differentiation_matrix(n) 
D = (u,u₀) -> M*(u .- u₀)/Δt
u₀ = 1.0; U = [u₀]; T = [0.0]; α = 0.9
for i = 0:N-1
  u = nlsolve(u->F(u,u₀,α),u₀*ones(n)).zero
  u₀ = u[end]
  append!(T,(i.+t)*Δt)
  append!(U,u)
end
plot(T, U, marker=:circle, legend=false)

A = [0    0    0    0    0; 
     1/3  0    0    0    0; 
     1/6  1/6  0    0    0;
     1/8  0    3/8  0    0;
     1/2  0   -3/2  2    0];
b = [1/6  0    0    2/3  1/6];

using LinearAlgebra, Plots
function rk_stability_plot(A,b)
  E = ones(length(b),1)
  r(λk) = abs.(1 .+ λk * b*((I - λk*A)\E))
  x,y = LinRange(-4,4,100),LinRange(-4,4,100)
  s = reshape(vcat(r.(x'.+im*y)...),(100,100))
  contour(x,y,s,levels=[1],legend=false)
end
rk_stability_plot(A,b)

i = 0:3
a = ((-i)'.^i.//factorial.(i))\[0;1;0;0]
b = ((-(i.+1))'.^i.//factorial.(i))\[1;0;0;0]
[a b]

function multistepcoefficients(m,n)
  s = length(m) + length(n) - 1
  A = (m.+1).^(0:s-1) .|> Rational
  B = (0:s-1).*((n.+1).^[0;0:s-2])
  c = -[A[:,2:end] B]\ones(Int64,s)
  [1;c[1:length(m)-1]], c[length(m):end] 
end

using Polynomials
function PECE(n,m)
  _,a = multistepcoefficients([0 1],hcat(1:n-1...))
  _,b = multistepcoefficients([0 1],hcat(0:n-1...))
  α(r) = a ⋅ r.^(1:n-1)/b[1]
  β(r) = b[2:end] ⋅ r.^(1:n-1)/b[1]
  z = [(c = [r-1; repeat([r + β(r)],m); α(r)];
    Polynomials.roots(Polynomial(c)))
    for r in exp.(im*LinRange(0,2π,200))]
  hcat(z/b[1]...)
end

z = vcat([PECE(2,i)[:] for i in 0:4]...)
scatter(z,label="",markersize=.5,aspect_ratio=:equal)

function pade(a,m,n)
  A = Rational(1)*Matrix(I(m+n+1))
  A[:,m+2:end] = [i>j ? -a[i-j] : 0 for i∈1:m+n+1, j∈1:n]
  pq = A\a[1:m+n+1]
  pq[1:m+1], [1; pq[m+2:end]]
end

m = 3; n = 2
a = [0; ((-1).^(0:m+n).//(1:m+n+1))]
(p,q) = pade(a,m,n)

S = n -> [(-1).^(i+j)*binomial(j,i) for i∈0:n, j∈0:n]
S(m)*p, S(n)*q

function SIR!(du,u,p,t)
  S,I,R = u
  β,γ = p
  du[1] = dS = -β*S*I
  du[2] = dI = +β*S*I - γ*I
  du[3] = dR = +γ*I
end

using DifferentialEquations
u₀ = [0.99; 0.01; 0]
tspan = (0.0,15.0)
p = (2.0,0.4)
problem = ODEProblem(SIR!,u₀,tspan,p)
solution = solve(problem)
plot(solution,labels=["susceptible" "infected" "recovered"])

using DifferentialEquations, Plots
function duffing!(dx,x,γ,t) 
  dx[1] = x[2]
  dx[2] = -γ*x[2]+x[1]-x[1]^3+0.3*cos(t)
end
problem = ODEProblem(duffing!,[1.0,0.0],(0.0,200.0),0.37)
solution = solve(problem,Vern7())
plot(solution, idxs=(1,2))

function solveIVP(f,u0,tspan)
  sol = solve(ODEProblem(f,u0,tspan))
  return(sol.u[end][1]) 
end

using DifferentialEquations, Roots
airy(y,p,x) = [y[2];x*y[1]]
domain = (-12.0,0.0);  bc = (1.0,1.0);  guess = 5
shoot_airy = (guess -> solveIVP(airy,[bc[1];guess],domain)-bc[2])
v = find_zero(shoot_airy, guess)

sol = solve(ODEProblem(airy,[bc[1],v],domain))
plot(sol)

m, n = [0 1 2 3 4], [1]
a, b = zeros(maximum(m)+1), zeros(maximum(n)+1) 
a[m.+1], b[n.+1] = multistepcoefficients(m,n)

λk = r -> (a ⋅ r.^-(1:length(a))) ./ (b ⋅ r.^-(1:length(b)))
r = LinRange(0.2,6,100)
plot([λk.(r) λk.(-r)],[r r],xlim=[-2,2])

Δx = 0.01; Δt = 0.01
ℓ = 1; x = -ℓ:Δx:ℓ; m = length(x)
Uⁿ = exp.(-8*x.^2); Uⁿ⁻¹ = Uⁿ  
ν = Δt/Δx^2; α = 0.5 + ν; γ = 0.5 - ν
B = ν*Tridiagonal(ones(m-1),zeros(m),ones(m-1))
B[1,2] *=2; B[end,end-1] *=2
@gif for i = 1:300
  global Uⁿ, Uⁿ⁻¹ = (B*Uⁿ+γ*Uⁿ⁻¹)/α, Uⁿ
  plot(x,Uⁿ,ylim=(0,1),label=:none,fill=(0, 0.3, :red))
end

function schroedinger(m,n,ε)
    x = LinRange(-4,4,m); Δx = x[2]-x[1]; Δt = 2π/n; V = x.^2/2
    ψ = exp.(-(x.-1).^2/2ε)/(π*ε)^(1/4)
    diags = 0.5im*ε*[1 -2 1]/Δx^2 .- im/ε*[0 1 0].*V
    D = Tridiagonal(diags[2:end,1], diags[:,2], diags[1:end-1,3])
    D[1,2] *= 2; D[end,end-1] *= 2
    A = I + 0.5Δt*D 
    B = I - 0.5Δt*D 
    for i ∈ 1:n 
      ψ = B\(A*ψ)
    end
    return ψ
end

ε = 0.3; m = 20000; eₓ=[]; eₜ=[]
N = floor.(Int,exp10.(LinRange(2,3.7,6)))
x = LinRange(-4,4,m)
ψₘ = -exp.(-(x.-1).^2/2ε)/(π*ε)^(1/4)
for n ∈ N
  x = LinRange(-4,4,n)
  ψₙ = -exp.(-(x.-1).^2/2ε)/(π*ε)^(1/4)
  append!(eₜ,norm(ψₘ - schroedinger(m,n,ε))/m)
  append!(eₓ,norm(ψₙ - schroedinger(n,m,ε))/n)
end
plot(2π./N,eₜ,shape=:circle, xaxis=:log, yaxis=:log)
plot!(8 ./N,eₓ,shape=:circle)

t = 0.5; n=100; m = 100
r = LinRange(0,2,m); Δr = r[2]-r[1]; Δt = t/n;
u₀ = @. tanh(32(1 - r)); u = u₀
d = @. [1 -2 1]/Δr^2  + (1/r)*[-1 0 1]/2Δr
D = Tridiagonal(d[2:end,1],d[:,2],d[1:end-1,3])
D[1,1:2] = [-4 4]/Δr^2; D[end,end-1:end] = [2 -2]/Δr^2 
A = I - 0.5Δt*D 
B = I + 0.5Δt*D  
for i = 1:n
 u = A\(B*u)
end
plot(r, u,label=:none,fill=(-1, 0.3, :red))

using Sundials
problem = ODEProblem((u,p,t)->D*u,u₀,(0,t))
method = CVODE_BDF(linear_solver=:Band,jac_upper=1,jac_lower=1)
solution = solve(problem,method)
plot(r, solution(t),label=:none,fill=(-1, 0.3, :red))

logitspace(x,n,p) = x*atanh.(LinRange(-p,p,n))/atanh(p)

function laplacian(x)
  Δx = diff(x); Δx₁ = Δx[1:end-1]; Δx₂ = Δx[2:end]
  d₋ = @.  2/[Δx₁*(Δx₁+Δx₂); Δx₁[1]^2]
  d₀ = @. -2/[Δx₂[end].^2; Δx₁*Δx₂; Δx₁[1].^2]  
  d₊ = @.  2/[Δx₂[end].^2; Δx₂*(Δx₁+Δx₂) ]
  Tridiagonal(d₋,d₀,d₊)
end

function heat_equation(x,t,n,u)
  Δt = t/n
  D² = laplacian(x)
  A = I - 0.5Δt*D²
  B = I + 0.5Δt*D²
  for i ∈ 1:n
    u = A\(B*u)
  end
  return u
end

ϕ = (x,t,s) -> exp.(-s*x.^2/(1+4*s*t))/sqrt(1+4*s*t)
m = 40; t = 15
x = LinRange(-20,20,m)
plot(x,ϕ(x,t,10),label="exact",width=3)
u₁ = heat_equation(x,t,m,ϕ(x,0,10))
plot!(x,u₁,shape=:circle,label="equal")
x = logitspace(20,m,0.999)
u₂ = heat_equation(x,t,m,ϕ(x,0,10))
plot!(x,u₂,shape=:circle,label="logit")

L = 16; m = 400; Δx = L/m
T = 4; n = 1600; Δt = T/n
x = LinRange(-L/2,L/2,m)'
u = @. tanh(x^4 - 16*(2*x^2-x'^2)) 
D = Tridiagonal(ones(m-1),-2ones(m),ones(m-1))/Δx^2
D[1,2] *= 2; D[end,end-1] *= 2
A = I + 0.5Δt*D
B = I - 0.5Δt*D 
f = (u,Δt) -> @. u/sqrt(u^2 - (u^2-1)*exp(-50*Δt))
u = f(u,Δt/2)
anim = Animation()
for i = 1:n
  (i%10)==1 && (plot(Gray.(u),border=:none); frame(anim))
  u = (B\(A*(B\(A*u))'))'
  (i<n) && (u = f(u,Δt))
end
u = f(u,Δt/2); Gray.(u)
gif(anim, "allencahn.gif", fps = 30)

U = (x,t) -> @. t<2 ? 
   (0<x<t)*(x/t) + (t<x<1+t/2) : (x/t)*(0<x<sqrt(2t))

m = 80; x = LinRange(-1,3,m); Δx = x[2]-x[1]; j = 1:m-1
n = 100; Lₜ = 4; Δt = Lₜ/n
f = u -> u.^2/2; df = u -> u
u = (x.>=0).*(x.<=1) 
anim = @animate for i = 1:n 
  α = 0.5*max.(abs.(df(u[j])),abs.(df(u[j.+1])))
  F = (f(u[j])+f(u[j.+1]))/2 - α.*(u[j.+1]-u[j])
  global u -=  Δt/Δx*[0;diff(F);0]
  plot(x,u, fill = (0, 0.3, :blue))
  plot!(x,U(x,i*Δt), fill = (0, 0.3, :red))  
  plot!(legend=:none, ylim=[0,1])
end
gif(anim, "burgers.gif", fps = 15)

δ = u -> diff(u,dims=1)
ϕ = t -> (abs(t)+t)./(1+abs(t))
fixnan(u) = isnan(u)||isinf(u) ? 0 : u
θ = δu -> fixnan.(δu[1:end-1,:]./δu[2:end,:])
∂ₓ(u) = (δu=δ(u);[[0 0];δu[2:end,:].*ϕ.(θ(δu));[0 0]])
F = u -> [u[:,1].*u[:,2] u[:,1].*u[:,2].^2+0.5u[:,1].^2]

m = 100; x = LinRange(-.5,.5,m); Δx = x[2]-x[1]
T = 0.4; n = ceil(T/(Δx/2)); Δt = (T/n)/2; 
U = [0.8*(x.<0).+0.2 zero(x)] 
Uⱼ = view(U,1:m-1,:); Uⱼ₊₁ = view(U,2:m,:)
anim = @animate  for i = 1:n
  U° = U-0.5*Δt/Δx*∂ₓ(F(U))
  Uⱼ₊₁ .= (Uⱼ+Uⱼ₊₁)/2 - δ(∂ₓ(U))/8 - Δt/Δx*δ(F(U°)) 
  U° = U-0.5*Δt/Δx*∂ₓ(F(U))
  Uⱼ .= (Uⱼ+Uⱼ₊₁)/2 - δ(∂ₓ(U))/8 - Δt/Δx*δ(F(U°))
  plot(x,U[:,1], fill = (0, 0.3, :blue))
  plot!(legend=:none, ylim=[0,1])
end
gif(anim, "dambreak.gif", fps = 15)

m=10; x=LinRange(0,1,m); h=x[2]-x[1]
α = fill(2/h-2h/3,m); α[[1,m]] /= 2; β = fill(-1/h-h/6,m-1)
A = SymTridiagonal(α,β)
b = [-2h^3/3; -4h^3/3 .- 8h*x[2:m-1].^2; -4h+8h^2/3-2h^3/3+1]
u = A\b
s = -16 .+ 8x.^2 .+ 15csc(1).*cos.(x)
plot(x,s,marker=:o,alpha=0.5); plot!(x,u,marker=:o,alpha=0.5)

m = 8; x = LinRange(0,1,m+2); h = x[2]-x[1]
σ = (a,b,c) -> Tridiagonal(fill(a,m-1),fill(b,m),fill(c,m-1))/h^3
M = [σ(-12,24,-12) σ(-6,0,6);σ(6,0,-6) σ(2,8,2)];
b = [384h*ones(m);zeros(m)]
u = M\b
s = 16*(x.^4 - 2x.^3 + x.^2)
plot(x,[s [0;u[1:m];0]],marker=:o,alpha=0.5)
s = 16*(x.^4 - 2x.^3 + x.^2)
plot(x,[s [0;u[1:m];0]],marker=:o,alpha=0.5)

using FFTW, DifferentialEquations
m = 128; x = (1:m)/m*(2π) .- π
ξ = im*[0:(m÷2); (-m÷2+1):-1]
f = (u,p,t) -> -real(ifft(ξ.*fft(0.5u.^2)))
u = solve(ODEProblem(f,exp.(-x.^2),(0.0,2.0)),DP5());

@gif for t=0:.02:2
  plot(x,u(t),ylim=(-0.2,1.2))
end

using FFTW, DifferentialEquations
ϕ = (x,x₀,c) -> 0.5c*sech(sqrt(c)/2*(x-x₀))^2
L = 30; T = 2.0; m = 256
x = (1:m)/m*L .- L/2
u₀ = ϕ.(x,-4,4) + ϕ.(x,-9,9)
iξ = im*[0:(m÷2);(-m÷2+1):-1]*2π/L
F = plan_fft(u₀)

G = t -> exp.(-iξ.^3*t)
f = (w,_,t) -> -G(t) .\ (3iξ .* (F * (F \ (G(t).*w) ).^2) )

w = solve(ODEProblem(f,F*u₀,(0,T)),DP5())
u = t -> real(F\(w(t).*G(t)))

@gif for t=0:.01:2
  plot(x,u(t),ylim=(0,5),legend=:none)
end

using FFTW, Images
ϵ = 1; m = 256; ℓ = 100; n = 2000; Δt=100/n
U = (rand(m,m).>.5) .- 0.5
ξ = [0:(m÷2);(-m÷2+1):-1]*2π/ℓ
D²= -ξ.^2 .- (ξ.^2)'
E = exp.(-(D².+1).^2*Δt)
F = plan_fft(U)
f = U -> U./sqrt.(U.^2/ϵ + exp(-Δt*ϵ)*(1 .- U.^2/ϵ))
for i=1:600
   U = f(F\(E.*(F*f(U))))
end
Gray.(clamp01.((real(U).+1)/2))

using FFTW, DifferentialEquations, SciMLOperators
ℓ = 8; T = 8.0; m = 256; ε² = 0.01; α = 1
ξ = [0:(m÷2);(-m÷2+1):-1]*2π/ℓ
Δ = -(ξ.^2 .+ ξ'.^2); Δ² = Δ.^2
u₀ = randn(m,m)
F = plan_fft(u₀)
L = DiagonalOperator((-ε²*Δ²+α*Δ)[:])
N = (u,_,_) -> (v = F\reshape(u,m,m); Δ.*(F*(v.^3-(1+α)*v)))[:]
problem = SplitODEProblem(L,N,(F*u₀)[:],(0,T))
solution = solve(problem,ETDRK4(),dt = 0.1)
u = t -> real(F\reshape(solution(t),m,m))

using Interact, ColorSchemes, Images
@manipulate for t in slider(0:0.1:8; value=0, label="time")
  get(colorschemes[:binary], u(t), :extrema)  
end

using FFTW, DifferentialEquations, SciMLOperators
ℓ = 50; T = 200.0; n = 128; u₀ = randn(n,n)
x = LinRange(-ℓ/2,ℓ/2,n+1)[2:end]
ξ = [0:(n÷2);(-n÷2+1):-1]*2π/ℓ
Δ = -(ξ.^2 .+ ξ'.^2); Δ² = Δ.^2
F = plan_fft(u₀)
L = DiagonalOperator((-Δ-Δ²)[:])
N = (u,_,_) -> 
  ( v = reshape(u,n,n); 
    v = -0.5*abs.((F\(im*ξ.*v)).^2 + (F\(im*ξ'.*v)).^2);
    w = (F*v)[:]; w[1] = 0.0;  return w )
problem = SplitODEProblem(L,N,(F*u₀)[:],(0,T))
solution = solve(problem, ETDRK4(), dt=0.05, saveat=0.5);
u = t -> real(F\reshape(solution(t),n,n))

using Interact, ColorSchemes, Images
@manipulate for t in slider(0:0.5:T; value=0, label="time")
  imresize(get(colorschemes[:magma], -u(t), :extrema),(256,256))  
end