library(plot3D)
library(pracma)
library(rmatio)
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=""))
}

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

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

mat = read.mat("nt_toolbox/data/ml-prostate.mat")
A = mat$A
class_names = mat$class_names[[1]]

head(A)

A = A[sample(dim(A)[1]),]

X = A[,1:(dim(A)[2] - 2)]
y = A[,dim(A)[2] - 1]
c = A[,dim(A)[2]]

n = dim(X)[1]
p = dim(X)[2]
print(n)
print(p)

I0 = (c==1) # train
I1 = (c==0) # test
n0 = sum(I0)
n1 = n - n0
X0 = X[I0,]
y0 = y[I0]
X1 = X[I1,]
y1 = y[I1]

mX0 = apply(X0, 2, mean)
sX0 = apply(X0, 2, std)
X0 = sweep(X0, 2, mX0,"-")
X0 = sweep(X0, 2, sX0,"/")
X1 = sweep(X1, 2, mX0,"-")
X1 = sweep(X1, 2, sX0,"/")

m0 = mean(y0)
y0 = y0 - m0
y1 = y1 - m0

options(repr.plot.width=4, repr.plot.height=4)
C = t(X0) %*% X0
image(1:p, 1:p, C, xlab="", ylab="", col=topo.colors(5), ylim=c(p + 0.5, 0.5))

options(repr.plot.width=4, repr.plot.height=4)
barplot(t(t(X0) %*% y0), xlab='', ylab='', col="darkblue")

svd_decomp = svd(X0)
U = svd_decomp$u
s = svd_decomp$d
V = svd_decomp$v
X0r = X0 %*% 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)

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(X0, col="blue", diag.panel=panel.hist)

options(repr.plot.width=4, repr.plot.height=4)
plot(X0r[,1], X0r[,2], col="blue", pch=16, xlab="", ylab="")

options(repr.plot.width=5, repr.plot.height=3)
for (i in 1:p)
{
    plot(X0r[,i], y0, col=i+1, xlab="", ylab="", pch=16)
}

w = solve(t(X0) %*% X0) %*% t(X0) %*% y0

options(repr.plot.width=4, repr.plot.height=4)
plot( X1 %*% w, type="l", col="blue", ylim=c(min(y1), max(y1)), xlab="", ylab="")
lines(y1,  type="l", col="orange")
legend("topright", legend=c("X1*w", "y1"), col=c("blue", "orange"), pch="-")

E = sqrt(sum((X1 %*% w-y1)^2 ) / n1)
print(paste('Relative prediction error:', E))

lambda = .2 * norm(X0)**2
w = solve(t(X0) %*% X0 + lambda * base::diag(p)) %*% t(X0) %*% y0
w1 = t(X0) %*% solve(X0 %*% t(X0) + lambda * base::diag(n0)) %*% y0
print(paste('Error (should be 0):', round(norm(w-w1)/norm(w), 3)))

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

## Insert your code here.

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

## Insert your code here.

J = function(w,Lambda){1/2 * norm(X0 %*% w - y0)**2 + Lambda * base::norm(as.matrix(w),"1")}

Soft = function(x,s){pmax(abs(x) - s, 0) * sign(x)}

options(repr.plot.width=4, repr.plot.height=4)
t = seq(-5,5,length=201)
plot(t, Soft(t,2), xlab="", ylab="", col=4, type="l") 

tau = 1.5 / base::norm(X0,"2")**2

lmax = max(abs(t(X0) %*% y0))
Lambda = lmax /10

w = rep(0, p)

C = t(X0) %*% X0
u = t(X0) %*% y0
ISTA = function(w,Lambda,tau){Soft(w - tau * (C %*% w - u), Lambda*tau)}
w = ISTA(w,Lambda,tau)

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

## Insert your code here.

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

## Insert your code here.

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

## Insert your code here.

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