options(warn=-1) # turns off warnings, to turn on: "options(warn=0)"
library(plot3D)
library(pracma)
library(grid)

# Importing the libraries
for (f in list.files(path="nt_toolbox/toolbox_general/", pattern="*.R")) {
    source(paste("nt_toolbox/toolbox_general/", f, sep=""))
}
for (f in list.files(path="nt_toolbox/toolbox_signal/", pattern="*.R")) {
    source(paste("nt_toolbox/toolbox_signal/", f, sep=""))
}

Xm = function(X){as.matrix(X - rep(colMeans(X), rep.int(nrow(X), ncol(X))))}
Cov = function(X){data.matrix(1. / (n - 1) * t(Xm(X)) %*% Xm(X))}

n = 1000 # number of sample
p = 2 # dimensionality
omega = 1.5 * 2.5 #offset
n1 = n/2
X = rbind(randn(n1,2), randn(n1,2) + rep(1, n1) * omega)
y = c(rep(1, n1), rep(-1, n1))

options(repr.plot.width=5, repr.plot.height=5)

for (i in c(-1, 1))
{
    I = (y==i)
    plot(X[I,1], X[I,2], col=(i + 3), xlim=c(min(X[,1]), max(X[,1])),  
         ylim=c(min(X[,2]), max(X[,2])), xlab="", ylab="", pch=16)
    par(new=TRUE)
}

cols = c(2,4)
legend("topright", legend=c(-1, 1), col=cols, pch=16)

options(repr.plot.width=7, repr.plot.height=6)
t = seq(-3, 3, length=255)
plot(t, log(1 + exp(t)), type="l", col=2)
#plot(t, t > 0)
lines(t, t > 0, col=3)
lines(t, pmax(t, 0), col=4)
legend("topleft", legend=c('Binary', 'Logistic', 'Hinge'), col=c(3,2,4), pch="-")

L = function(s,y){1/n * sum( log(1 + exp(-s * y)))}
E = function(w,X,y){L(X %*% w, y)}

theta = function(v){1 / (1 + exp(-v))}
nablaL = function(s, r){ - 1/n * y * theta(-s * y)}
nablaE = function(w,X,y){t(X) %*% nablaL(X %*% w,y)}

AddBias = function(X){cbind(X, rep(1, dim(X)[1]))}

w = rep(0, p + 1)
dim(w) = c(p+1, 1)

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

tau_max = 2/(1/4 * base::norm(AddBias(X), "2")**2)
print(tau_max)

source("nt_solutions/ml_3_classification/exo1.R")

## Insert your code here.

q = 201
tx = seq(min(X[,1]), max(X[,1]), length=q)
ty = seq(min(X[,2]),max(X[,2]), length=q)
B = as.vector(meshgrid(ty, tx)$X)
A = as.vector(meshgrid(ty, tx)$Y)
G = matrix(c(A, B), nrow=length(A), ncol=2)

Theta = theta(AddBias(G) %*% w)
dim(Theta) = c(q, q)

color = function(x){rev(cm.colors(x))}
image(tx,ty, Theta, xlab="", ylab="", col=color(10), xaxt="n", yaxt="n")
      par(new=TRUE)
for (i in c(-1, 1))
{
    I = (y==i)
    plot(X[I,1], X[I,2], col=(i + 3), xlim=c(min(X[,1]), max(X[,1])),  
         ylim=c(min(X[,2]), max(X[,2])), xlab="", ylab="", pch=16, xaxt="n", yaxt="n")
    par(new=TRUE)
}

cols = c(2,4)
legend("topright", legend=c(-1, 1), col=cols, pch=16)

source("nt_solutions/ml_3_classification/exo2.R")

## Insert your code here.

source("nt_solutions/ml_3_classification/exo3.R")

## Insert your code here.

distmat = function(X,Z)
{
    dist1 = diag(X %*% t(X))
    dist2 = diag(Z %*% t(Z))
    n1 = dim(X)[1]
    n2 = dim(Z)[1]
    out = matrix(0, n1, n2)
    for (i in 1:n1)
    {
        for (j in 1:n2)
        {
            out[i,j] = dist1[i] + dist2[j]
        }
    }
    out = out - 2 * X %*% t(Z)

    return(out)
}

kappa = function(X,Z,sigma){exp( -distmat(X,Z)/(2*sigma^2))}

n = 1000
p = 2
t = 2 * pi * rand(n/2,1)
R = 2.5
r = R * (1 + .2 * rand(n/2,1)) # radius
X1 = cbind(cos(t) * r, sin(t) * r)
X = rbind(randn(n/2, 2), X1)
y = c(rep(1, n/2), rep(-1, n/2))

options(repr.plot.width=5, repr.plot.height=5)

for (i in c(-1, 1))
{
    I = (y==i)
    plot(X[I,1], X[I,2], col=(i + 3), xlim=c(min(X[,1]), max(X[,1])),  
         ylim=c(min(X[,2]), max(X[,2])), xlab="", ylab="", pch=16)
    par(new=TRUE)
}

cols = c(2,4)
legend("topright", legend=c(-1, 1), col=cols, pch=16)

sigma = 1
K = kappa(X, X, sigma)
image(K, col=color(10), ylim=c(1, 0))

F = function(h,K,y){L(K %*% h, y)}
nablaF = function(h,K,y){t(K) %*% nablaL(K %*% h,y)}

source("nt_solutions/ml_3_classification/exo4.R")

## Insert your code here.

q = 201
tmax = 3.5
t = seq(-tmax, tmax, length=q)
B = as.vector(meshgrid(t)$X)
A = as.vector(meshgrid(t)$Y)
G = matrix(c(A, B), nrow=length(A), ncol=2)
Theta = theta(kappa(G,X,sigma) %*% h)
dim(Theta) = c(q, q)

image(t,t, Theta, xlab="", ylab="", col=color(10), xaxt="n", yaxt="n")
par(new=TRUE)
for (i in c(-1, 1))
{
    I = (y==i)
    plot(X[I,1], X[I,2], col=(i + 3), xlim=c(min(X[,1]), max(X[,1])),  
         ylim=c(min(X[,2]), max(X[,2])), xlab="", ylab="", pch=16, xaxt="n", yaxt="n")
    par(new=TRUE)
}

cols = c(2,4)
legend("topright", legend=c(-1, 1), col=cols, pch=16)

source("nt_solutions/ml_3_classification/exo5.R")

## Insert your code here.

source("nt_solutions/ml_3_classification/exo6.R")

## Insert your code here.

LSE0 = function(S){log(apply(exp(S), 1, sum))}

max2 = function(S){apply(S, 1, max)}
LSE = function(S){LSE0(S - max2(S)) + max2(S)}

# check equality of LSE and LSE0
S = randn(4,5)
norm( LSE(S)-LSE0(S))

SM0 = function(S){exp(S) / apply(exp(S), 1, sum)}

SM = function(S){SM0(S-max2(S))}

# Check equality of SM and SM0
norm( SM(S)-SM0(S) )

file_name = 'nt_toolbox/data/digits.csv'
A = read.table(file_name, sep=",", head=FALSE)
#A = A[sample(dim(A)[1]),]
X = as.matrix(A[,1:(dim(A)[2] - 1)])
y = A[,dim(A)[2]]

n = dim(X)[1]
p = dim(X)[2]
CL = sort(unique(y)) # list of classes.
k = length(CL)

options(repr.plot.width=6, repr.plot.height=6)
par(mar = rep(2, 4))
par(mfrow=c(k, 5))

q = 5
for (i in 1:k)
{
    I = which(y==CL[i])
    for (j in 1:q)
    {
        
        f = as.numeric(X[I[j],])
        f = f / max(f)
        dim(f) = c(sqrt(p), sqrt(p))
        image(-f[,ncol(f):1], col=gray(c(0,1)), xaxt="n", yaxt="n", bty="n")
    }
}

svd_decomp = svd(Xm(X))
U = svd_decomp$u
s = svd_decomp$d
V = svd_decomp$v
X0r = Xm(X) %*% V

options(repr.plot.width=4, repr.plot.height=4)
plot(s, type="o", col=4, ylab="", xlab="", pch=16)

options(repr.plot.width=6, repr.plot.height=6)

plot_multiclasses(X, y, 2)

plot_multiclasses(X,y, 3)

D = matrix(0, nrow=n, ncol=k)
for (i in 1:n)
{
    for (j in 1:k)
    {
        if (y[i] == (j - 1)){D[i, j] = 1}
    }
}

E = function(W){(1./n) * (sum(LSE(X %*% W)) - (c(X %*% W) %*% c(D)))}

nablaE = function(W){1./n * t(X) %*% ( SM(X %*% W) -  D)}

source("nt_solutions/ml_3_classification/exo7.R")

## Insert your code here.

svd_decomp = svd(Xm(X))
U = svd_decomp$u
D = svd_decomp$d
V = t(svd_decomp$v)
Z = Xm(X) %*% t(V)
M = max(abs(c(Z)))
        
q = 201
t = seq(-M, M, length=q)
B = as.vector(meshgrid(t, t)$X)
A = as.vector(meshgrid(t, t)$Y)

G0 = matrix(0, nrow=q*q, ncol=2)
G0[,1:2] = c(B, A)
G = G0 %*% V[1:2,] + matrix(rep(apply(X, 2, mean), q*q), nrow=q*q, ncol=dim(X)[2], byrow=TRUE)

Theta = SM(G %*% W)
dim(Theta) = c(q, q, k)

par(mfrow=c(3, 4))

for (i in 1:k)
{
    image(t(Theta[,,i]), main=paste("Class", i), ylim=c(1, 0), xlim=c(0, 1), axes=FALSE)
}

plot_multiclasses_bis = function(X, Y, dim)
{
    
    svd_decomp = svd(Xm(X))
    U = svd_decomp$u
    D = svd_decomp$d
    V = svd_decomp$v
    Z = Xm(X) %*% V
    classes = sort(unique(y))
    nb_classes = length(classes)

    for (i in classes)
    {
        I = (y==i)
        color = col[,i + 1]
        plot(Z[I,1], Z[I,2], col=rgb(color[1], color[2], color[3]), xlim=c(min(Z[,1]), max(Z[,1])),  
             ylim=c(min(Z[,2]), max(Z[,2])), xlab="", ylab="", pch=16, axes=FALSE)
        par(new=TRUE)
    }
}

col = c(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)
col = matrix(col, nrow=3, ncol=k)

R = array(0, dim=c(q,q,3))

for (i in 1:k)
{
    for (a in 1:3)
    {
        R[,,a] = R[,,a] + adjust(Theta[,,i]) * col[a,i]
    }
}

options(repr.plot.width=5, repr.plot.height=5)
plot_multiclasses_bis(X, y, 2)
par(new=TRUE)
grid.raster(aperm(R, c(1, 2, 3)), interpolate=FALSE)
plot_multiclasses_bis(X, y, 2)

source("nt_solutions/ml_3_classification/exo8.R")

## Insert your code here.