compressedsensing.py

import cvxpy as cp
import numpy as np
import matplotlib.pyplot as plt
import random
from tqdm import tqdm


def main():
  tree_list = []
  with open("wonderland-tree.txt", "r") as f:
    for line in f:
      tree_list += [int(c) for c in line.strip()]

  x = np.array(tree_list)
  n = x.size
  print("Sparsity: %f" % (x.sum() / n))

  """
  A = np.random.normal(size=(n, n))

  b = A[:700] @ x
  x_rec = recover(b, A[:700])
  print("x_700 is close enough to x: %r" % np.allclose(x, x_rec))

  r_min = bin_search(A, x, n) # 614
  print("r_min: %d" % r_min)

  norms = []
  for i in tqdm(range(-10, 3)):
    x_rec = recover(A[:r_min + i] @ x, A[:r_min + i])
    norms.append(np.linalg.norm(x - x_rec, ord=1))

  print(norms)

  fig, ax = plt.subplots()
  ax.set_title("L1-norm versus Compression Index")
  ax.set_ylabel("L1-norm of Difference between recovered x and x")
  ax.set_xlabel("Compressed Dimension relative to r*")
  ax.plot(range(-10, 3), norms)
  plt.show()

  rps = []
  for p in np.arange(0.011, 0.021, 0.001):
    print("\n\np: ", p)
    rs = []

    for i in range(5):
      A = np.zeros((n, n))

      for i in range(n):
        for j in range(n):
          prob = random.random()

          if prob <= p / 2:
            A[i, j] = - 1
          elif prob <= p:
            A[i, j] = 1

      r_min = bin_search(A, x, n)
      rs.append(r_min)

    r_hat = sum(rs) / len(rs) # 628.0
    print(r_hat)
    rps.append(r_hat)

  print("rps: ", rps)
  fig, ax = plt.subplots()
  ax.set_title("Average Minimum r for Probability p")
  ax.set_ylabel("Minimum r")
  ax.set_xlabel("p")
  ax.plot(rps)
  plt.show()

  """


def recover(b, A):
  x = cp.Variable(1200)
  constraints = [A @ x == b]
  obj = cp.Minimize(cp.atoms.norm1(x))
  prob = cp.Problem(obj, constraints)
  prob.solve()

  return x.value


def bin_search(A, x, n):
  l, r = 0, 1200
  while l + 1 < r:
    m = int((r + l) / 2)
    x_rec = recover(A[:m] @ x, A[:m])

    if np.linalg.norm(x - x_rec, ord=1) < 0.001:
      r = m
    else:
      l = m

  return r

if __name__ == "__main__":
  main()