using PyPlot
using NtToolBox

eta = 8;

f = x -> (x[1]^2 + eta*x[2]^2) / 2;

f([2,3])

include("NtToolBox/src/ndgrid.jl")
t = linspace(-.7,.7,101)
(u, v) = meshgrid(t,t)
F = (u.^2 + eta.*v.^2) ./ 2
contourf(t, t, F, 35);

function Grad_f(x)
    return[x[1]; eta*x[2]]
end;

Grad_f([1,2]);

tau = 1.8/max(eta,1); tau

nbiter = 10
x = [1, 1]  # initial estimate of the solution
for iter in 1:nbiter  # iter goes from 1 to nbiter
    x = x - tau.*Grad_f(x)
end
x   # to display x, like in Matlab. Use print(x) if this is not the last command of the cell, else nothing is displayed

f(x)

function GradDescent(Grad_f, x0, nbiter, tau)
    x = x0
    for iter in 1:nbiter  # iter goes from 1 to nbiter
        x = x - tau.*Grad_f(x)   # x has type 'array'. Like in C, one can also write x-=tau*Grad_f(x)
    end
    return x
end;

GradDescent(Grad_f, [1, 1], 10, tau)

function GradDescentArray(Grad_f, x0, nbiter, tau)
    x = x0
    xlist = [x0[1] x0[2]]  
    for iter in 1:nbiter
        x = x - tau*Grad_f(x)     
        xlist = [xlist; [x[1] x[2]]]
    end
    return xlist
end;

xarray = GradDescentArray(Grad_f,[0.6,0.6],10,0.225)

xarray[2, 1]

xarray[2, :]

xarray[:, 1]

mapslices(f, xarray, 2)

plot(0:size(xarray)[1] - 1, log10(mapslices(f, xarray, 2)), "o-");

contourf(t, t, F, 35)
xarray = GradDescentArray(Grad_f, [0.6, 0.6], 30, 0.225)
plot(xarray[:, 1], xarray[:, 2], "w.-")

name = "NtToolBox/src/data/lena.png"
xsharp = load_image(name)*256
println(size(xsharp))
println("The size of the image is $(size(xsharp,1)) x $(size(xsharp,2)).")
println("The range of the pixel values is [$(minimum(xsharp)), $(maximum(xsharp))].") # the range of the pixel values is [0.0, 1.0] because load_image scales the image.

xsharp = readdlm("xsharp");

imageplot(xsharp)
set_cmap("gray")
colorbar()       # displays the color bar close to the image
subplots_adjust(top = 0.75)
title("This is Lena");

function D(x)
    vdiff = [diff(x,1) ; zeros(1,size(x,2))] # concatenates along the rows
    hdiff = [diff(x,2) zeros(size(x,1),1)]   # concatenates along the columns
    return cat(3, vdiff, hdiff) # combination along a third dimension
end;

v = D(xsharp);

fig = figure(figsize=(16,7))
suptitle("The two images of vertical and horizontal finite differences")
subplot(1,2,1)
imageplot(v[:, :, 1])
title("Image of vertical finite differences")
set_cmap("gray")
colorbar()

subplot(1,2,2)
imageplot(v[:, :, 2])
colorbar()
title("Image of horizontal finite differences")
set_cmap("gray")

imshow(sqrt(sum(v.^2, 3)[:, :]))
colorbar()      
subplots_adjust(top = 0.75);

Dadj = v -> [-v[1,:,1]'; -diff(v[1:end-1,:,1],1); v[end-1,:,1]'] + 
            [-v[:,1,2] -diff(v[:,1:end-1,2],2)  v[:,end-1,2]];

sum((D(xsharp).^2)) - sum(Dadj(D(xsharp)).*xsharp)

(N1,N2) = size(xsharp)
noiselevel = 20/N1
y = xsharp + noiselevel * randn(N1,N2);

figure(figsize=(11,11))
imshow(y, interpolation="nearest")
set_cmap("gray")
colorbar()       
subplots_adjust(top=0.75)
title("This is a noisy version of Lena");

function Grad_f(x, y, Lambda)
    x - y + Lambda .* Dadj(D(x))
end;

include("NtSolutions/optim_1bis_gradient_descent/exo1.jl");

function GradDescent(Grad_f, y, Lambda, x0, nbiter, tau)
    # put your code here.
    x = x0;
    for iter in 1 : nbiter  
        x = x - tau.*Grad_f(x, y, Lambda)
    end
end;

include("NtSolutions/optim_1bis_gradient_descent/exo2.jl");

Lambda = 8
xdenoised = GradDescent(Grad_f, y, Lambda, y, 100, 1.9/(1+8*Lambda));

figure(figsize=(10,10))
imshow(xdenoised)
set_cmap("gray")
title("Denoised image");

mask = rand(size(xsharp,1),size(xsharp,2)).>0.9;

fig= figure(figsize=(16,7)) # one figure with two horizontal subfigures
subplot(1,2,1)
imshow(mask)
set_cmap("gray")
title("The binary mask")
subplot(1,2,2)
imshow(mask.*xsharp)
set_cmap("gray")
title("The available pixel values are displayed, the missing pixels are in black");

y = mask.*xsharp;

include("NtSolutions/optim_1bis_gradient_descent/exo3.jl");

function Proj_Omega(x, y, mask)
    x[mask] = y[mask] # parameters of functions are given by reference, so the content of x can be modified
    return x
end;

include("NtSolutions/optim_1bis_gradient_descent/exo4.jl");

tau = 1.9 / 8
nbiter = 300;

xrestored = ProjGradDescent(Grad_f, Proj_Omega, y, mask, y, nbiter, tau);

figure(figsize=(10,10))
imshow(xrestored)
set_cmap("gray")
title("Denoised image");
<script>
  $(document).ready(function(){
      $('div.prompt').hide();
  });
</script>