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

# 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=""))
}

data(iris)
A = iris
# Transforming the species to numerical categories :
A$Species = as.numeric(A$Species)

A = A[sample(nrow(A)),]

X = as.matrix(A[,1:(dim(A)[2] - 1)])
y = as.matrix(A[,dim(A)[2]])
y = y - min(y) + 1

n = dim(X)[1]
p = dim(X)[2]
k = max(y)

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

panel.hist = function(x, ...)
{
    usr = par("usr"); on.exit(par(usr))
    par(usr = c(usr[1:2], 0, 1.5) )
    h = hist(x, plot = FALSE)
    breaks = h$breaks; nB <- length(breaks)
    y = h$counts
    y = y/max(y)
    rect(breaks[-nB], 0, breaks[-1], y, col = "darkblue")
}
pairs(X, col=c("red","blue","green")[y], diag.panel=panel.hist)

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

image(c(1:4), c(1:4), Cov(X), xlab='', ylab='')

svd_decomp = svd(Xm(X))
U = svd_decomp$u
D = svd_decomp$d
V = svd_decomp$v

Z = Xm(X) %*% V

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

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

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

cols = c(2:(k + 1))
legend("topright", legend=c(1:k), col=cols, pch=16)

for (i in 1:k)
{
    I = (y==i)
    if (i == 1)
    {
        scatter3D(Z[I,1], Z[I,2], Z[I,3], col=i + 1, xlim=c(min(Z[,1]), max(Z[,1])),  
         ylim=c(min(Z[,2]), max(Z[,2])), zlim=c(min(Z[,3]), max(Z[,3])), xlab="", ylab="",
          zlab="", pch=16)
    }
    else
    {
         scatter3D(Z[I,1], Z[I,2], Z[I,3], col=i + 1, xlim=c(min(Z[,1]), max(Z[,1])),  
         ylim=c(min(Z[,2]), max(Z[,2])), zlim=c(min(Z[,3]), max(Z[,3])), xlab="", ylab="",
          zlab="", pch=16, add=TRUE)
    }
}

cols = c(2:(k + 1))
legend("topright", legend=c(1:k), col=cols, pch=16)

n0 = round(.5 * n)
n1 = n - n0
X0 = X[1:n0,]
y0 = y[1:n0]
X1 = X[(n0+1):dim(X)[1],]
y1 = y[(n0+1):dim(y)[1]]

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)
}

i = 1 
x = X1[i,] # could be any point
x = t(as.matrix(x, c(length(x))))
D = distmat(X0,x)

ys = y0[order(D)]

# Defining an hist function
hist = function(y, classes)
{
    if (length(dim(y)) == 2)
    {
        out = matrix(0, dim(y)[2], length(classes))
        for (i in 1:dim(y)[2])
        {
            out[i,] = hist(y[,i], classes)
        }
        return(out)
    }
    else
    {
        count = c()
        for (i in classes)
        {
            count = c(count, sum(y == i))
        }
        return(count)
    }
}

R = 5
h = hist(ys[1:R], 1:k) / R
c = which.max(h)
print(paste('c(x) =',c,"[true class =", y1[i], "]"))

options(repr.plot.width=3, repr.plot.height=3)
Rlist = round(c(0.05, 0.1, 0.5, 1) * n0)

for (R in Rlist)
{
    h = hist(ys[1:R], 1:k) / R
    barplot(h, main=paste("Value of R :", R), col=4)
}

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

## Insert your code here.

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

## Insert your code here.

I = sample(n)
I = I[1:k]
C = X[I,]

D = distmat(X, C)
yb = apply(D, 1, which.min)

for (i in 1:k)
{
    I = yb==i
    plot(Z[I,1], Z[I,2], col=(i+1), xlim=c(min(Z[,1]), max(Z[,1])),  
         ylim=c(min(Z[,2]), max(Z[,2])), xaxt='n', yaxt='n', 
         ylab="", xlab="", bty="n")
    par(new=TRUE)
}

CV = (C - rep(colMeans(X), rep.int(nrow(C), ncol(X)))) %*% V

for (i in 1:k)
{
    plot(CV[i, 1], CV[i, 2], , xlim=c(min(Z[,1]), max(Z[,1])),  
         ylim=c(min(Z[,2]), max(Z[,2])), xaxt='n', yaxt='n', 
         ylab="", xlab="", lwd=5, bty="n")
    par(new=TRUE)
}

for (l in 1:k)
{
    C[l,] = apply(X[yb==l,], 2, mean)
}

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

## Insert your code here.

# Insert your code here.

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

## Insert your code here.