using PyPlot
using NtToolBox
# using Autoreload
# arequire("NtToolBox")

remove_diag = C -> C - diagm(diag(C))
Correlation = Phi -> remove_diag(abs(Phi'*Phi));

mu = Phi -> maximum(Correlation(Phi));

Coh = (Phi, k) -> (k*mu(Phi))/(1 - (k - 1)*mu(Phi));

normalize = Phi -> Phi ./ repeat(sqrt(sum(Phi.^2, 1)), inner = [size(Phi)[1], 1]);


PhiRand = (P, N) -> normalize(randn(P, N))
Phi = PhiRand(250, 1000);

c = mu(Phi)
println(@sprintf("Coherence:    %.2f", c))
println(@sprintf("Sparsity max: %d", floor(1/2*(1 + 1/c))))

include("NtSolutions/sparsity_6_l1_recovery/exo1.jl")

## Insert your code here.

supp   = s -> find(abs(s) .> 1e-5)
cosupp = s -> find(abs(s) .< 1e-5);

PsiI = (Phi,I) -> Phi[:, setdiff(1:size(Phi)[2], I)]' * pinv(Phi[:,I])';

F = (Phi, s) -> norm(PsiI(Phi, supp(s))*s[supp(s)], Inf);

erc =(Phi, I) -> norm(PsiI(Phi,I), Inf);

g = (C, I) -> sum(C[:,I], 2)
werc_g = (g, I, J) -> maximum(g[J]) / (1 - maximum(g[I]))
werc = (Phi, I) -> werc_g(g(Correlation(Phi), I), I, setdiff(1:size(Phi)[2], I) );

include("NtSolutions/sparsity_6_l1_recovery/exo2.jl")

## Insert your code here.

N = 2000
P = N - 10
Phi = PhiRand(N, P)
s = zeros(N)
s[1:6] = 1
I = supp(s)
k = length(I)

println(@sprintf("N = %d, P = %d, |I| = %d", N, P, k))
println(@sprintf("F(s)     = %.2f",  F(Phi, s)))
println(@sprintf("ERC(I)   = %.2f",  erc(Phi, I)))
println(@sprintf("w-ERC(s) = %.2f", werc(Phi, I)))
println(@sprintf("Coh(|s|) = %.2f", Coh(Phi, k)))

include("NtSolutions/sparsity_6_l1_recovery/exo3.jl")

## Insert your code here.

N = 600
P = Int(N/2)
Phi = PhiRand(P, N)
klist = Array{Int64,1}(round(linspace(1, P/7., 20)))
ntrials = 60
proba = zeros(length(klist), 3)

for i in 1:length(klist)
    proba[i, 1:3] = 0
    k = Int(klist[i])
    for j in 1:ntrials
        s = zeros(N)
        I = randperm(N)
        I = I[1:k]
        l = randn(k, 1)
        s[I] = l./abs(l)
        proba[i, 1] = proba[i, 1] + (F(Phi, s) .< 1)
        proba[i, 2] = proba[i, 2] + (erc(Phi, I) .< 1)
        proba[i, 3] = proba[i, 3] + (werc(Phi, I) .> 0).*(werc(Phi, I) .< 1)
    end
end
        
figure(figsize = (8, 5))
plot(klist, proba/ntrials, linewidth = 2)
xlabel("k")
legend(["F <1", "ERC <1", "w-ERC <1"])
title(@sprintf("N = %d, P = %d", N, P))
show()

minmax = v -> (1 - minimum(v), maximum(v) - 1)
ric = A -> minmax(eig(A'*A)[1]);

include("NtSolutions/sparsity_6_l1_recovery/exo4.jl")

## Insert your code here.

include("NtSolutions/sparsity_6_l1_recovery/exo5.jl")

## Insert your code here.

sigma = 6
g = x -> (1 - x.^2./sigma^2).*exp(-x.^2./(2*sigma^2));

include("NtToolBox/src/ndgrid.jl")
P = 1024
(Y, X) = meshgrid(0:P-1, 0:P-1)
Phi = normalize(g((X - Y + P/2.)%P-P/2.));

eta = 2
N = P/eta
Phi = Phi[:, 1:eta:end];

c = Phi'*Phi
c = abs(c[:, Int(size(c)[2]/2)])


figure(figsize = (8, 5))
plot(c[Base.div(length(c), 2) - 50 : Base.div(length(c), 2) + 50], ".-")
show()

twosparse = d -> circshift([1; zeros(d, 1); -1; zeros(Int(N) - d - 2, 1)], round(N/2 - d/2));

x0 = twosparse(50)


figure(figsize = (10, 7))
subplot(2, 1, 1)
plot(x0[:], "r", linewidth = 2)
subplot(2, 1, 2)
plot((Phi*x0)[:], "b", linewidth = 2)
show()

include("NtSolutions/sparsity_6_l1_recovery/exo6.jl")

## Insert your code here.
<script>
  $(document).ready(function(){
      $('div.prompt').hide();
  });
</script>