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="")) } remove_diag = function(C){C - base::diag(base::diag(C))} Correlation = function(Phi){remove_diag(abs(t(Phi) %*% Phi))} mu = function(Phi){max(Correlation(Phi))} Coh = function(Phi, k){(k * mu(Phi)) / (1 - (k-1) * mu(Phi))} normalize = function(Phi){t(t(Phi) / sqrt(apply(Phi**2,2, sum)))} PhiRand = function(P, N){normalize(matrix( rnorm(P*N,mean=0,sd=1), P, N))} Phi = PhiRand(250, 1000) c = mu(Phi) print(paste("Coherence:", round(c, 2))) print(paste("Sparsity max:", floor(1/2*(1 + 1/c)))) source("nt_solutions/sparsity_6_l1_recovery/exo1.R") ##Insert your code here. where = function(I) { n = length(I) out = c() for (i in 1:n) { if (I[i] == TRUE) { out = c(out, i) } } return(out) } supp = function(s){where(abs(s) > 1e-5)} cosupp = function(s){where(abs(s) < 1e-5)} PsiI = function(Phi,I) { keep_indices = c() for (i in 1:dim(Phi)[2]) { if (!(i %in% I)) { keep_indices = c(keep_indices, i) } } return (t(Phi[,keep_indices]) %*% t(pinv(as.matrix(Phi[,I])))) } F = function(Phi,s){base::norm(PsiI(Phi, supp(s)) %*% s[supp(s)], type="I")} erc = function(Phi, I){base::norm(PsiI(Phi, I), type="I")} setdiff = function(n, I) {keep_indices = c() for (i in 1:n) { if (!(i %in% I)) { keep_indices = c(keep_indices, i) } } return(keep_indices) } g = function(C,I){apply(as.matrix(C[,I]), 1, sum)} werc_g = function(g,I,J){max(g[J])/(1 - max(g[I]))} werc = function(Phi,I){werc_g(g(Correlation(Phi), I), I, setdiff(dim(Phi)[2], I))} source("nt_solutions/sparsity_6_l1_recovery/exo2.R") ## Insert your code here. source("nt_solutions/sparsity_6_l1_recovery/exo3.R") #Insert your code here. minmax = function(v){c(1 - min(v), max(v) - 1)} ric = function(A){minmax(eigen(t(A) %*% A)$values)} source("nt_solutions/sparsity_6_l1_recovery/exo4.R") ## Insert your code here. source("nt_solutions/sparsity_6_l1_recovery/exo5.R") ## Insert your code here. sigma = 6 g = function(x){(1 - (x**2 / sigma**2)) * exp(-x**2/(2*sigma**2))} P = 1024 t = meshgrid(0:(P - 1),0:(P - 1)) Y = t$X X = t$Y Phi = normalize(g((X - Y + P/2.) %% P - P/2.)) eta = 2. N = P/eta Phi = Phi[,seq(1, dim(Phi)[2], by=eta)] c = t(Phi) %*% Phi c = abs(c[,dim(c)[2]/2]) options(repr.plot.width=5, repr.plot.height=5) plot(c[(length(c)/2 - 50):((length(c)/2) + 50)], type="l", xlab="", ylab="", col=4) roll <- function( x , n ){ if( n == 0 ) return( x ) c( tail(x,n) , head(x,-n) ) } twosparse = function(d){roll(c(1, rep(0, d), -1, rep(0, N - d- 2)), as.integer(N/2 - d/2))} x0 = twosparse(50) plot(x0, type="l", col=2, ylab="", xlab="") plot(Phi %*% x0, type="l", col=4, ylab="", xlab="") source("nt_solutions/sparsity_6_l1_recovery/exo6.R") #Insert your code here.