set.seed(1234)
library(ADMM)
delta <- 0.1
n <- 50; p <- 200;

probs_0 <- NULL;
probs_1 <- NULL;
probs_4 <- NULL;

for(k in seq(1, 25, by=2)) {

  prob_0 <- 0;
  prob_1 <- 0;
  prob_4 <- 0;

  for(j in 1:100) {
    A <- matrix(rnorm(n*p, 0, 1), nrow=n, ncol=p)
    xstark <- c(rnorm(k, 0, 1), rep(0, 200-k))
    b <- A %*% xstark

    x_0 <- admm_bp(A, b)$fit()$beta

    D <- diag(1/(abs(x_0)+delta)[,1])
    gamma_bar <- admm_bp(A %*% solve(D), b)$fit()$beta
    x_1 <- solve(D) %*% gamma_bar

    D <- diag(1/(abs(x_1)+delta)[,1])
    gamma_bar <- admm_bp(A %*% solve(D), b)$fit()$beta
    x_2 <- solve(D) %*% gamma_bar

    D <- diag(1/(abs(x_2)+delta)[,1])
    gamma_bar <- admm_bp(A %*% solve(D), b)$fit()$beta
    x_3 <- solve(D) %*% gamma_bar

    D <- diag(1/(abs(x_3)+delta)[,1])
    gamma_bar <- admm_bp(A %*% solve(D), b)$fit()$beta
    x_4 <- solve(D) %*% gamma_bar

    if(max(abs(x_0 - xstark)) < 0.01) prob_0 = prob_0 + 1
    if(max(abs(x_1 - xstark)) < 0.01) prob_1 = prob_1 + 1
    if(max(abs(x_4 - xstark)) < 0.01) prob_4 = prob_4 + 1
  }

  probs_0 <- c(probs_0, prob_0/100)
  probs_1 <- c(probs_1, prob_1/100)
  probs_4 <- c(probs_4, prob_4/100)
}