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

library(pracma)
library(imager)

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

n = 128
f0 = load_image("nt_toolbox/data/lena.png")
f0 = t(f0[,])
f0 = rescale(f0[(256 - n / 2 + 1):(256 + n / 2), (256 - n / 2 + 1) : (256 + n / 2)])

options(repr.plot.width=4, repr.plot.height=4)
imageplot(f0, 'Image f_0')

rho = 0.7

set.seed(1)
Omega = c(matrix(0, n, n))
sel = sample(n**2)
Omega[sel[1:as.integer(rho * n**2)]] = 1
Omega = matrix(Omega, n, n)

Phi = function(f, Omega){f * (1 - Omega)}

y = Phi(f0, Omega)

imageplot(y, 'Observations y')

SoftThresh = function(x, T)
{
    if (!is.null(dim(x)))
    {
        return(x * pmax(1 - T / abs(x), 1e-10 * array(1, dim=dim(x)), array(0, dim=dim(x))))
    }
    else
    {
        return(x * pmax(1 - T / abs(x), 1e-10 * rep(1, length(x)), rep(0, length(x))))
    }
}

x = seq(-1, 1, length=1000)

plot(x, SoftThresh(x, .5), col=4, type="l", ylab="", xlab="")

Jmax = log2(n) - 1
Jmin = (Jmax - 3)

Psi = function(a){perform_wavelet_transf(a, Jmin, -1, ti=0)}
PsiS = function(f){perform_wavelet_transf(f, Jmin, +1, ti=0)}

SoftThreshPsi = function(f, T){Psi(SoftThresh(PsiS(f), T))}

imageplot(clamp(SoftThreshPsi(f0, 0.1)))

lambd = .03

ProjC = function(f, Omega){Omega * f + (1 - Omega) * y}

fSpars = y

fSpars = ProjC(fSpars, Omega)

fSpars = SoftThreshPsi(fSpars, lambd)

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

# Insert your code here.

# Insert your code here.

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

# Insert your code here.

J = Jmax - Jmin + 1
u = c(4**(-J), 4**(floor(-1 * seq(J + 2./3,1,by=-1 /3.) + 1)))
U = array(0, dim=c(n, n, length(u)))
for (i in 1:length(u))
{
    U[,,i] = u[i]
}
U = aperm(U, c(3, 1, 2))

lambd = .01

Xi = function(a){perform_wavelet_transf(a, Jmin, -1, ti=1)}
PsiS = function(f){perform_wavelet_transf(f, Jmin, + 1, ti=1)}
Psi = function(a){Xi(a/U)}

tau = 1.9 * min(u)

a = U * PsiS(fSpars)

fTI = Psi(a)
a = a + tau * PsiS(Phi(y - Phi(fTI, Omega), Omega))

a = SoftThresh(a, lambd * tau)

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

## Insert your code here.

## Insert your code here.

## Insert your code here.

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

## Insert your code here.

HardThresh = function(x, t){x * (abs(x) > t)}

x = seq(-1, 1, length=1000)

plot(x, HardThresh(x, .5), col=4, type="l", xlab="", ylab="")

niter = 500

lambda_list = seq(1, 0, length=niter)

fHard = y

fHard = ProjC(fHard, Omega)

fHard = Xi(HardThresh(PsiS(fHard), tau * lambda_list[1]))

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

## Insert your code here.