using PyPlot
using NtToolBox

#SetAR = @(ar)set(gca, 'PlotBoxAspectRatio', [1 ar 1], 'FontSize', 10);
Xm = X -> X - repeat(mean(X,1), outer=(size(X,1), 1))
Cov = X -> Xm(X)'*Xm(X)
dotp = (u,v) ->sum(u[:].*v[:]);

n = 1000 # number of sample
p = 2 # dimensionality
omega = [1 .5]*5 # offset 
X = [randn(div(n,2),2); randn(div(n,2),2)+ones(div(n,2),1)*omega]
y = [ones(div(n,2),1);-ones(div(n,2),1)];

plot_multiclasses(X,y,disp_dim=2,ms=5);

t = linspace(-3,3,255);
plot(t, log(1+exp(t)), c="red",lw=2, label="Logistic")
plot(t, t.>0, c="orange",lw=2, label="Hinge")
plot(t,max(t,0), c="blue", label="Binary")
axis("tight"); legend(loc="best");

L = (s,y) -> 1/n * sum( log( 1 + exp(-s.*y) ) );
E = (w,X,y) -> L(X*w,y);

theta = v -> 1 ./ (1+exp(-v));
nablaL = (s,r) -> - 1/n * y.* theta(-s.*y);
nablaE = (w,X,y) -> X'*nablaL(X*w,y);

AddBias = X -> [X ones(size(X,1),1)];

w = zeros(p+1,1);

tau = .8; # here we are using a fixed tau
w = w - tau * nablaE(w,AddBias(X),y);

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

# Insert your code here.

q = 201; tx = linspace(minimum(X[:,1]),maximum(X[:,1]),q); ty = linspace(minimum(X[:,2]),maximum(X[:,2]),q)
B,A = meshgrid( ty,tx )
G = [A[:] B[:]];

Theta = theta(AddBias(G)*w);
Theta = reshape(Theta, q, q);

imshow(Theta'[:,end:-1:1], extent=[minimum(tx), maximum(tx), minimum(ty), maximum(ty)]);
plot_multiclasses(X,y);

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

# Insert your code here.

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

# Insert your code here.

distmat = (X,Z) -> broadcast(+,sum(X'.*X',1)',sum(Z'.*Z',1))-2*(X*Z');

kappa = (X,Z,sigma) -> exp( -distmat(X,Z)/(2*sigma^2) );

n = 1000; p = 2
t = 2*pi*rand(div(n,2),1)
R = 2.5
r = R*(1 + .2*rand(div(n,2),1)); # radius
X1 = [cos(t).*r sin(t).*r]
X = [randn(div(n,2),2); X1]
y = [ones(div(n,2),1);-ones(div(n,2),1)];

plot_multiclasses(X,y,disp_dim=2;disp_legend=1);
axis("off");

sigma = 1
K = kappa(X,X,sigma);

F = (h,K,y) -> L(K*h,y);
nablaF = (h,K,y) -> K'*nablaL(K*h,y);

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

# Insert your code here.

q = 201
tmax = 3.5
t = linspace(-tmax,tmax,q)
B,A = meshgrid( t,t )
G = [A[:] B[:]]
Theta = reshape( theta(kappa(G,X,sigma)*h) , q,q);

imshow(Theta'[:,end:-1:1], extent=[-tmax, tmax, -tmax, tmax], cmap = get_cmap("jet"));
plot_multiclasses(X,y,disp_legend = 0);
clim(0,1)
#caxis([0 1]);
axis("off");

include("NtSolutions/ml_3_classification/exo5.jl");

# Insert your code here.

include("NtSolutions/ml_3_classification/exo6.jl");

# Insert your code here.

LSE = S -> log( sum(exp(S), 2) );

max2 = S -> repeat(mapslices(maximum, S,2), outer=(1, size(S,2)))
LSE_stab = S -> LSE( S-max2(S) ) + mapslices(maximum, S,2);

SM = S -> exp(S) ./ repeat( sum(exp(S),2), outer=(1, size(S,2)));

SM_stab = S -> SM(S-max2(S));

A = readdlm("NtToolBox/src/data/digits.csv", ',')
A = A[randperm(size(A,1)),:]
X = A[:,1:end-1]; y = A[:,end];

n,p = size(X);
CL = unique(y); # list of classes.
k = length(CL);

q = 5
for i=1:k
    I = find(y.==CL[i]);
    for j=1:q
        f = reshape(X[I[j],:], Int(sqrt(p)), Int(sqrt(p)))';
        subplot(q,k, (j-1)*k+i );
        imshow(-f, extent=[0,1,0,1], cmap=get_cmap("gray")); axis("image"); axis("off")
    end
end


plot_multiclasses(X,y,disp_dim = 2, ms = 5, disp_legend = 1);

plot_multiclasses(X,y,disp_dim=3);

D = Int.(repeat(CL[:]', outer=(n,1)) .== repeat(y, outer=(1,k))); #double

E = W -> 1/n*( sum(LSE(X*W)) - dotp(X*W,D)  );

nablaE = W -> 1/n * X'* ( SM_stab(X*W) -  D );

include("NtSolutions/ml_3_classification/exo7.jl");

W = zeros(p,k)
Elist = []
tau = .01
niter = 500
for i=1:niter
    W = W - tau * nablaE(W)
    append!(Elist, E(W))
end
plot(1:niter, Elist, lw=2)
axis("tight");

# Insert your code here.

U,D,V = svd(Xm(X),thin=true)
Z = Xm(X) * V
M = maximum(abs(Z[:]))
q = 201
t = linspace(-M,M,q)
B,A = meshgrid(t,t)
G = zeros(q*q,p)
G[:,1:2] = [A[:]; B[:]]
G = G*V' + repeat(mean(X,1), outer=(q*q, 1));

Theta = SM(G*W)
Theta = reshape(Theta, q, q, k);

for i=1:k
    subplot(3,4,i)
    imshow(Theta[:,:,i]'[:,end:-1:1], extent = [0,1,0,1], cmap=get_cmap("jet"));
    title(string("Class ", i))
    axis("image"); axis("off")
end
tight_layout()

col = [ [1 0 0]; [0 1 0]; [0 0 1]; [0 0 0]; [0 1 1]; [1 0 1]; [1 1 0];
    [1 .5 .5]; [.5 1 .5]; [.5 .5 1]  ]'
R = zeros(q,q,3)
for i=1:k
    for a=1:3
        R[:,:,a] = R[:,:,a] + Theta[:,:,i] .* col[a,i]
    end
end

imshow(permutedims(R, [2, 1, 3]), extent = [-tmax, tmax, -tmax, tmax]); 
plot_multiclasses3(X,y,disp_dim=2,ms=3,disp_legend=1);
axis("off");

include("NtSolutions/ml_3_classification/exo8.jl");

# Insert your code here.