using PyPlot, Interact, LinearAlgebra
arc(args...; kws...) = gca().add_patch.(matplotlib.patches.Arc(args...; kws...))

A = [1   1
    -1 1/4]
fig = figure()

Θ = range(0,2π, length=200)
V = [[cos(θ),sin(θ)] for θ in Θ]
U = [A*v for v in V]
Vx, Vy = first.(V), last.(V)
Ux, Uy = first.(U), last.(U)

# @manipulate for θ in slider(0:1:90, value=29)
for θ in 9:10:49
    display(
    withfig(fig) do
        ϕ = deg2rad(θ) # convert to radians
        v₁, v₂ = [cos(ϕ), sin(ϕ)], [-sin(ϕ), cos(ϕ)] # orthonormal inputs
        u₁, u₂ = A*v₁, A*v₂ # outputs
        
        # arrows and labels for v₁, v₂, u₁, u₂
        for (v,s,c) in ((v₁,"v₁","red"), (v₂,"v₂","red"), (u₁,"Av₁","blue"), (u₂,"Av₂","blue"))
            arrow(0,0,v..., color=c, width=0.04, length_includes_head=true)
            text(((norm(v)+0.1)*normalize(v))..., s, color=c)
        end
        
        # show the angle between u₁ and u₂
        uangle = rad2deg(acos((u₁⋅u₂)/(norm(u₁)*norm(u₂))))
        arc([0,0], 0.5,0.5, angle=rad2deg(atan(u₁[2],u₁[1])), theta1=0, theta2=uangle, color="blue")
        p = 0.4*normalize(u₁+u₂)
        text(p..., "$(round(uangle,sigdigits=2))°", color="blue")
        if round(uangle,sigdigits=2) == 90
            arc(p + [0.2,0.05], 0.5, 0.23, color="blue")
        end
        
        # plot dashed lines for all possible v₁, v₂, u₁, u₂
        
        plot(Vx,Vy, "r--")
        plot(Ux,Uy, "b--")
        
        axis("square")
        xlim(-2,2)
        ylim(-2,2)
    end
    )
end

U, σ, V = svd(A)
σ # the ellipsoid semi-axes

V # columns of V are the orthogonal "input" basis v₁, v₂

for (v,s,c) in ((V[:,1],"v₁","red"),  (V[:,2],"v₂","red"),
                (U[:,1],"u₁","blue"), (U[:,2],"u₂","blue"))
    arrow(0,0,v..., color=c, width=0.04, length_includes_head=true)
    text(((norm(v)+0.15)*normalize(v))..., s, color=c)
end
axis("square")
xlim(-1.5,1.5)
ylim(-1.5,1.5)
title(L"the singular vectors of $A$")

rad2deg(acos(-V[1,1]))  # rotation angle of the basis, in degrees

using Images, FileIO

picture = download("http://web.mit.edu/jfrench/Public/gstrang.png")
pimage = load(picture)

p = Float64.(channelview(pimage)) # convert to an array
summary(p)

pr,pg,pb = p[1,:,:],p[2,:,:],p[3,:,:]
Ur,σr,Vr = svd(pr)
Ug,σg,Vg = svd(pg)
Ub,σb,Vb = svd(pb);

subplot(1,3,1)
imshow(cat(pr,0*pg,0*pb, dims=3))
subplot(1,3,2)
imshow(cat(0*pr,pg,0*pb, dims=3))
subplot(1,3,3)
imshow(cat(0*pr,0*pg,pb, dims=3))

semilogy(σr, "r.")
semilogy(σg, "g.")
semilogy(σb, "b.")
title("exponential decay of singular values of Strang image")
ylabel(L"\sigma_k")
xlabel(L"k")
legend(["red","green","blue"])

# clip x to [0,1] so that imshow doesn't complain about rounding errors
# that lead to values slightly outside this range.
clip01(x) = max(min(x, 1), 0)

fig = figure()
# @manipulate for k=slider(1:40,value=1)
for k in 1:5:26
    display(
    withfig(fig) do
        title("rank $k approximation")
        p̂r = clip01.(Ur[:,1:k]*Diagonal(σr[1:k])*Vr[:,1:k]')
        p̂g = clip01.(Ug[:,1:k]*Diagonal(σg[1:k])*Vg[:,1:k]')
        p̂b = clip01.(Ub[:,1:k]*Diagonal(σb[1:k])*Vb[:,1:k]')
        imshow(cat(p̂r,p̂g,p̂b, dims=3))
    end
    )
end

pics = Dict{String,Array}()
pics["Ukraine"] = load(download("https://upload.wikimedia.org/wikipedia/commons/thumb/4/49/Flag_of_Ukraine.svg/640px-Flag_of_Ukraine.svg.png"))
pics["USA"] = load(download("https://upload.wikimedia.org/wikipedia/en/thumb/a/a4/Flag_of_the_United_States.svg/640px-Flag_of_the_United_States.svg.png"))
# pics["Klingon"] = load(download("https://images-na.ssl-images-amazon.com/images/I/51ibu5dAb9L._SY550_.jpg"))
display(pics["Ukraine"])
display(pics["USA"])
# display(pics["Klingon"])

fig = figure()
# @manipulate for flagimg in collect(keys(pics)), k=slider(1:40,value=1)
for flagimg in collect(keys(pics)), k=1:10
  display(
    withfig(fig) do
        title("rank-$k approximation of $flagimg flag")
        p = float.(channelview(pics[flagimg])) # convert to an array
        pr,pg,pb = p[1,:,:],p[2,:,:],p[3,:,:]
        Ur,σr,Vr = svd(pr)
        Ug,σg,Vg = svd(pg)
        Ub,σb,Vb = svd(pb)
        p̂r = clip01.(Ur[:,1:k]*Diagonal(σr[1:k])*Vr[:,1:k]')
        p̂g = clip01.(Ug[:,1:k]*Diagonal(σg[1:k])*Vg[:,1:k]')
        p̂b = clip01.(Ub[:,1:k]*Diagonal(σb[1:k])*Vb[:,1:k]')
        imshow(cat(p̂r,p̂g,p̂b, dims=3))
    end
  )
end

B = [1 2
     1 2.01]
σ = svdvals(B) # returns just the σ values

σ[1]/σ[2]

cond(B) # same thing: the condition number of B

Θ = range(0,2π, length=200)
V = [[cos(θ),sin(θ)] for θ in Θ]
U = [B*v for v in V]
Vx, Vy = first.(V), last.(V)
Ux, Uy = first.(U), last.(U)
plot(Vx,Vy, "r--")
plot(Ux,Uy, "b--")

axis("square")
xlim(-4,4)
ylim(-4,4)
legend(["input circle", "output ellipse"])
title("inputs and outputs of ill-conditioned B")

U,σ,V = svd(B)
B̃ = σ[1] * U[:,1] * V[:,1]' # "best" rank-1 approximation to B

B̃ - B

x = B \ [1, 2]

B \ [1, 1]

B \ [1, 1.01]