Compressed Sensing Image

get_ipython().ast_node_interactivity = 'all'
import os
import matplotlib.pyplot as plt
import numpy as np
import matplotlib
import math
from PIL import Image
import scipy.fftpack
import scipy.optimize
import random
from sklearn.linear_model import Lasso
matplotlib.rcParams['figure.dpi'] = 150

def show(img):
    _ = plt.imshow(img, cmap='gray')
    #_ = plt.colorbar()

def img_open(path):
    img = Image.open(path)
    img = img.convert('L')
    return np.array(img, dtype=np.uint8)

target = img_open("/home/leo/Downloads/test.jpg")

W = target.shape[1]
H = target.shape[0]

show(target)

received = int((W * H) * 0.05)
not_received = (W * H) - received

k = [True] * received + [False] * not_received
random.shuffle(k)
random.shuffle(k)
random.shuffle(k)
k = np.array(k)

b = target.T.flat[k].astype(float)

show(target * k.reshape(target.shape))

def dct2(x):
    return scipy.fftpack.dct(scipy.fftpack.dct(x.T, norm='ortho', axis=0).T, norm='ortho', axis=0)

def idct2(x):
    return scipy.fftpack.idct(scipy.fftpack.idct(x.T, norm='ortho', axis=0).T, norm='ortho', axis=0)

def evaluate(x, g, step):
    """An in-memory evaluation callback."""

    # we want to return two things: 
    # (1) the norm squared of the residuals, sum((Ax-b).^2), and
    # (2) the gradient 2*A'(Ax-b)

    # expand x columns-first
    x2 = x.reshape((W, H)).T

    # Ax is just the inverse 2D dct of x2
    Ax2 = idct2(x2)

    # stack columns and extract samples
    Ax = Ax2.T.flat[k].reshape(b.shape)

    # calculate the residual Ax-b and its 2-norm squared
    Axb = Ax - b
    fx = np.sum(np.power(Axb, 2))

    # project residual vector (k x 1) onto blank image (ny x nx)
    Axb2 = np.zeros(x2.shape)
    Axb2.T.flat[k] = Axb # fill columns-first

    # A'(Ax-b) is just the 2D dct of Axb2
    AtAxb2 = 2 * dct2(Axb2)
    AtAxb = AtAxb2.T.reshape(x.shape) # stack columns

    # copy over the gradient vector
    np.copyto(g, AtAxb)

    return fx

from pylbfgs import owlqn

Xat2 = owlqn(W*H, evaluate, None, 5)


# transform the output back into the spatial domain
Xat = Xat2.reshape(W, H).T # stack columns
Xa = idct2(Xat)

show(Xa)

Iteration 1:
  fx = 1052846684.137251, xnorm = 1306.498417, gnorm = 64961.181239, step = 0.000015, k = 1, ls = 1

Iteration 2:
  fx = 378771374.740408, xnorm = 20718.937800, gnorm = 30439.029848, step = 1.000000, k = 2, ls = 1

Iteration 3:
  fx = 141900150.440163, xnorm = 40912.278714, gnorm = 3168.067844, step = 1.000000, k = 3, ls = 1

Iteration 4:
  fx = 137338582.827604, xnorm = 40297.907389, gnorm = 2663.509726, step = 1.000000, k = 4, ls = 1

Iteration 5:
  fx = 124123051.811049, xnorm = 39657.024794, gnorm = 4937.634041, step = 1.000000, k = 5, ls = 1

Iteration 6:
  fx = 117342007.507820, xnorm = 40076.228997, gnorm = 3376.995410, step = 1.000000, k = 6, ls = 1

Iteration 7:
  fx = 115136991.186309, xnorm = 48271.568151, gnorm = 8218.653274, step = 1.000000, k = 7, ls = 1

Iteration 8:
  fx = 110287628.397871, xnorm = 39629.126026, gnorm = 4814.367083, step = 1.000000, k = 8, ls = 1

Iteration 9:
  fx = 106708120.820831, xnorm = 47271.422723, gnorm = 6259.988157, step = 0.500000, k = 9, ls = 2

Iteration 10:
  fx = 100011073.232687, xnorm = 42174.376915, gnorm = 1850.074230, step = 0.500000, k = 10, ls = 2

Iteration 11:
  fx = 98521861.551988, xnorm = 47184.519170, gnorm = 4395.402566, step = 0.500000, k = 11, ls = 2

Iteration 12:
  fx = 95311661.101823, xnorm = 43586.628175, gnorm = 1361.503705, step = 1.000000, k = 12, ls = 1

Iteration 13:
  fx = 92091089.253206, xnorm = 47876.134240, gnorm = 3024.506824, step = 1.000000, k = 13, ls = 1

Iteration 14:
  fx = 88993950.366124, xnorm = 47611.130730, gnorm = 1845.838158, step = 1.000000, k = 14, ls = 1

Iteration 15:
  fx = 80679789.880600, xnorm = 51225.865527, gnorm = 1806.838550, step = 1.000000, k = 15, ls = 1

Iteration 16:
  fx = 80430189.070877, xnorm = 65552.480300, gnorm = 5643.934467, step = 1.000000, k = 16, ls = 1

Iteration 17:
  fx = 76479621.166977, xnorm = 55075.133108, gnorm = 2784.487806, step = 1.000000, k = 17, ls = 1

Iteration 18:
  fx = 70442781.512268, xnorm = 65036.793494, gnorm = 3037.624694, step = 0.250000, k = 18, ls = 3

Iteration 19:
  fx = 66308577.681624, xnorm = 64163.204384, gnorm = 1735.585194, step = 1.000000, k = 19, ls = 1

Iteration 20:
  fx = 64888293.389100, xnorm = 77773.066392, gnorm = 4287.672962, step = 0.500000, k = 20, ls = 2

Iteration 21:
  fx = 60926881.746555, xnorm = 70238.089090, gnorm = 1193.263490, step = 0.500000, k = 21, ls = 2

Iteration 22:
  fx = 56467129.590133, xnorm = 83010.213764, gnorm = 2950.675210, step = 1.000000, k = 22, ls = 1

Iteration 23:
  fx = 51384418.891695, xnorm = 82323.130568, gnorm = 1183.951568, step = 1.000000, k = 23, ls = 1

Iteration 24:
  fx = 48654496.810289, xnorm = 106734.319122, gnorm = 4341.616082, step = 1.000000, k = 24, ls = 1

Iteration 25:
  fx = 42719247.056787, xnorm = 95897.074118, gnorm = 1372.886405, step = 0.500000, k = 25, ls = 2

Iteration 26:
  fx = 40596869.463378, xnorm = 109036.673268, gnorm = 2757.796683, step = 0.500000, k = 26, ls = 2

Iteration 27:
  fx = 38041214.804251, xnorm = 101391.964330, gnorm = 1064.199873, step = 1.000000, k = 27, ls = 1

Iteration 28:
  fx = 36505941.937765, xnorm = 112647.199963, gnorm = 2244.358723, step = 0.250000, k = 28, ls = 3

Iteration 29:
  fx = 32639426.691367, xnorm = 112401.969964, gnorm = 1045.079168, step = 1.000000, k = 29, ls = 1

Iteration 30:
  fx = 30713005.437708, xnorm = 122905.419279, gnorm = 1736.104673, step = 0.500000, k = 30, ls = 2

Iteration 31:
  fx = 27730875.878713, xnorm = 122861.734136, gnorm = 1044.490693, step = 1.000000, k = 31, ls = 1

Iteration 32:
  fx = 26737465.188275, xnorm = 128781.509862, gnorm = 1227.193544, step = 0.250000, k = 32, ls = 3

Iteration 33:
  fx = 24967720.048042, xnorm = 128997.698605, gnorm = 1201.675812, step = 1.000000, k = 33, ls = 1

Iteration 34:
  fx = 24409187.244297, xnorm = 134995.518871, gnorm = 1165.097512, step = 0.250000, k = 34, ls = 3

Iteration 35:
  fx = 23971800.216224, xnorm = 130276.893341, gnorm = 853.354207, step = 0.500000, k = 35, ls = 2

Iteration 36:
  fx = 23519620.841449, xnorm = 135321.573008, gnorm = 927.212359, step = 0.250000, k = 36, ls = 3

Iteration 37:
  fx = 22750107.508386, xnorm = 134047.263174, gnorm = 612.031080, step = 1.000000, k = 37, ls = 1

Iteration 38:
  fx = 22012616.716144, xnorm = 138605.493301, gnorm = 751.082366, step = 1.000000, k = 38, ls = 1

Iteration 39:
  fx = 21590127.849869, xnorm = 136549.870394, gnorm = 549.016573, step = 1.000000, k = 39, ls = 1

Iteration 40:
  fx = 21247698.796309, xnorm = 139224.113592, gnorm = 553.108922, step = 1.000000, k = 40, ls = 1

Iteration 41:
  fx = 21004360.792569, xnorm = 137108.535380, gnorm = 420.531754, step = 1.000000, k = 41, ls = 1

Iteration 42:
  fx = 20802816.355241, xnorm = 139069.840718, gnorm = 416.734625, step = 0.500000, k = 42, ls = 2

Iteration 43:
  fx = 20570918.856432, xnorm = 138621.277167, gnorm = 315.438806, step = 1.000000, k = 43, ls = 1

Iteration 44:
  fx = 20429154.304982, xnorm = 140785.294031, gnorm = 456.719704, step = 1.000000, k = 44, ls = 1

Iteration 45:
  fx = 20290566.707010, xnorm = 139351.107952, gnorm = 281.262887, step = 0.500000, k = 45, ls = 2

Iteration 46:
  fx = 20196559.475057, xnorm = 140383.452146, gnorm = 292.117666, step = 1.000000, k = 46, ls = 1

Iteration 47:
  fx = 20117869.588810, xnorm = 139608.578336, gnorm = 219.476842, step = 1.000000, k = 47, ls = 1

Iteration 48:
  fx = 20020985.220258, xnorm = 140171.507354, gnorm = 212.246325, step = 1.000000, k = 48, ls = 1

Iteration 49:
  fx = 19951841.742712, xnorm = 140009.169427, gnorm = 168.290580, step = 1.000000, k = 49, ls = 1

Iteration 50:
  fx = 19890090.199447, xnorm = 140555.943020, gnorm = 188.386554, step = 1.000000, k = 50, ls = 1

Iteration 51:
  fx = 19867931.784446, xnorm = 140144.184583, gnorm = 142.823364, step = 1.000000, k = 51, ls = 1

Iteration 52:
  fx = 19833323.431126, xnorm = 140558.721908, gnorm = 139.424281, step = 1.000000, k = 52, ls = 1

Iteration 53:
  fx = 19801859.273860, xnorm = 140379.338724, gnorm = 109.053803, step = 1.000000, k = 53, ls = 1

Iteration 54:
  fx = 19775891.617246, xnorm = 140646.406696, gnorm = 130.030237, step = 1.000000, k = 54, ls = 1

Iteration 55:
  fx = 19765967.063161, xnorm = 140331.503621, gnorm = 102.435312, step = 1.000000, k = 55, ls = 1

Iteration 56:
  fx = 19765653.678841, xnorm = 140816.597265, gnorm = 161.495685, step = 1.000000, k = 56, ls = 1

Iteration 57:
  fx = 19743234.987051, xnorm = 140473.510267, gnorm = 82.223188, step = 0.500000, k = 57, ls = 2

Iteration 58:
  fx = 19730216.386376, xnorm = 140555.466947, gnorm = 76.317794, step = 1.000000, k = 58, ls = 1

Iteration 59:
  fx = 19716418.914154, xnorm = 140620.438055, gnorm = 69.233841, step = 1.000000, k = 59, ls = 1

Iteration 60:
  fx = 19705499.245321, xnorm = 140777.307697, gnorm = 68.000840, step = 1.000000, k = 60, ls = 1

Iteration 61:
  fx = 19697999.793255, xnorm = 140774.356214, gnorm = 58.866487, step = 1.000000, k = 61, ls = 1

Iteration 62:
  fx = 19689127.938768, xnorm = 140813.199814, gnorm = 59.192170, step = 1.000000, k = 62, ls = 1

Iteration 63:
  fx = 19685131.747294, xnorm = 140811.291026, gnorm = 51.095914, step = 1.000000, k = 63, ls = 1

Iteration 64:
  fx = 19678124.376384, xnorm = 140833.785885, gnorm = 53.546542, step = 1.000000, k = 64, ls = 1

Iteration 65:
  fx = 19675841.296515, xnorm = 140838.999364, gnorm = 45.055664, step = 1.000000, k = 65, ls = 1

Iteration 66:
  fx = 19671572.653538, xnorm = 140847.244516, gnorm = 36.829948, step = 1.000000, k = 66, ls = 1

Iteration 67:
  fx = 19666923.334664, xnorm = 140866.665107, gnorm = 42.134361, step = 1.000000, k = 67, ls = 1

Iteration 68:
  fx = 19665491.167786, xnorm = 140871.660941, gnorm = 35.825103, step = 1.000000, k = 68, ls = 1

Iteration 69:
  fx = 19662628.894564, xnorm = 140886.982073, gnorm = 30.250060, step = 1.000000, k = 69, ls = 1

Iteration 70:
  fx = 19659683.651725, xnorm = 140901.671638, gnorm = 29.139925, step = 1.000000, k = 70, ls = 1

Iteration 71:
  fx = 19657698.071720, xnorm = 140915.351270, gnorm = 25.951732, step = 1.000000, k = 71, ls = 1

Iteration 72:
  fx = 19655941.249919, xnorm = 140934.446724, gnorm = 26.375729, step = 1.000000, k = 72, ls = 1

Iteration 73:
  fx = 19654838.315670, xnorm = 140938.756121, gnorm = 23.278485, step = 1.000000, k = 73, ls = 1

Iteration 74:
  fx = 19653452.006347, xnorm = 140944.989450, gnorm = 23.188835, step = 1.000000, k = 74, ls = 1

Iteration 75:
  fx = 19652613.727277, xnorm = 140954.343300, gnorm = 20.257541, step = 1.000000, k = 75, ls = 1

Iteration 76:
  fx = 19651506.448323, xnorm = 140958.305472, gnorm = 19.597548, step = 1.000000, k = 76, ls = 1

Iteration 77:
  fx = 19650720.020801, xnorm = 140967.789868, gnorm = 17.388838, step = 1.000000, k = 77, ls = 1

Iteration 78:
  fx = 19649920.611807, xnorm = 140974.220955, gnorm = 16.668420, step = 1.000000, k = 78, ls = 1

Iteration 79:
  fx = 19649328.171528, xnorm = 140978.482309, gnorm = 15.387213, step = 1.000000, k = 79, ls = 1

Iteration 80:
  fx = 19648784.806142, xnorm = 140982.496732, gnorm = 14.562073, step = 1.000000, k = 80, ls = 1

Iteration 81:
  fx = 19648323.844189, xnorm = 140988.729557, gnorm = 13.832167, step = 1.000000, k = 81, ls = 1

Iteration 82:
  fx = 19647927.792271, xnorm = 140993.187955, gnorm = 12.989534, step = 1.000000, k = 82, ls = 1

Iteration 83:
  fx = 19647572.407780, xnorm = 140996.862796, gnorm = 12.653570, step = 1.000000, k = 83, ls = 1

Iteration 84:
  fx = 19647284.663178, xnorm = 141001.419321, gnorm = 12.053505, step = 1.000000, k = 84, ls = 1

Iteration 85:
  fx = 19647020.388558, xnorm = 141003.800974, gnorm = 11.231340, step = 1.000000, k = 85, ls = 1

Iteration 86:
  fx = 19646765.813671, xnorm = 141007.242356, gnorm = 10.391320, step = 1.000000, k = 86, ls = 1

Iteration 87:
  fx = 19646522.619968, xnorm = 141009.760657, gnorm = 9.910231, step = 1.000000, k = 87, ls = 1

Iteration 88:
  fx = 19646324.892210, xnorm = 141013.639747, gnorm = 9.511209, step = 1.000000, k = 88, ls = 1

Iteration 89:
  fx = 19646157.087521, xnorm = 141016.055956, gnorm = 9.099476, step = 1.000000, k = 89, ls = 1

Iteration 90:
  fx = 19646006.782719, xnorm = 141019.478997, gnorm = 8.476928, step = 1.000000, k = 90, ls = 1

Iteration 91:
  fx = 19645864.631023, xnorm = 141021.476950, gnorm = 7.985786, step = 1.000000, k = 91, ls = 1

Iteration 92:
  fx = 19645735.439932, xnorm = 141024.459353, gnorm = 7.514659, step = 1.000000, k = 92, ls = 1

Iteration 93:
  fx = 19645623.907647, xnorm = 141026.388817, gnorm = 7.312401, step = 1.000000, k = 93, ls = 1

Iteration 94:
  fx = 19645529.187373, xnorm = 141026.653346, gnorm = 6.969928, step = 1.000000, k = 94, ls = 1

Iteration 95:
  fx = 19645444.691133, xnorm = 141028.189823, gnorm = 6.586964, step = 1.000000, k = 95, ls = 1

Iteration 96:
  fx = 19645363.516632, xnorm = 141028.900411, gnorm = 6.046951, step = 1.000000, k = 96, ls = 1

Iteration 97:
  fx = 19645284.076586, xnorm = 141030.599818, gnorm = 5.831677, step = 1.000000, k = 97, ls = 1

Iteration 98:
  fx = 19645220.873989, xnorm = 141031.890259, gnorm = 5.785639, step = 1.000000, k = 98, ls = 1

Iteration 99:
  fx = 19645169.622285, xnorm = 141033.980885, gnorm = 5.642054, step = 1.000000, k = 99, ls = 1

Iteration 100:
  fx = 19645122.035141, xnorm = 141034.738286, gnorm = 5.109277, step = 1.000000, k = 100, ls = 1

Iteration 101:
  fx = 19645072.242173, xnorm = 141036.358699, gnorm = 4.567606, step = 1.000000, k = 101, ls = 1

Iteration 102:
  fx = 19645020.019335, xnorm = 141036.489967, gnorm = 4.519958, step = 1.000000, k = 102, ls = 1

Iteration 103:
  fx = 19644983.388942, xnorm = 141038.422760, gnorm = 4.498863, step = 1.000000, k = 103, ls = 1

Iteration 104:
  fx = 19644953.369891, xnorm = 141038.620752, gnorm = 4.400341, step = 1.000000, k = 104, ls = 1

Iteration 105:
  fx = 19644925.451267, xnorm = 141040.202358, gnorm = 4.047175, step = 1.000000, k = 105, ls = 1

Iteration 106:
  fx = 19644896.995925, xnorm = 141040.049114, gnorm = 3.623814, step = 1.000000, k = 106, ls = 1

Iteration 107:
  fx = 19644867.309876, xnorm = 141041.569613, gnorm = 3.566308, step = 1.000000, k = 107, ls = 1

Iteration 108:
  fx = 19644845.808363, xnorm = 141041.353911, gnorm = 3.752888, step = 1.000000, k = 108, ls = 1

Iteration 109:
  fx = 19644830.124102, xnorm = 141042.468314, gnorm = 3.799478, step = 1.000000, k = 109, ls = 1

Iteration 110:
  fx = 19644813.995621, xnorm = 141042.155191, gnorm = 3.352568, step = 1.000000, k = 110, ls = 1

Iteration 111:
  fx = 19644798.018471, xnorm = 141042.946785, gnorm = 2.882305, step = 1.000000, k = 111, ls = 1

Iteration 112:
  fx = 19644777.187703, xnorm = 141043.144443, gnorm = 2.618539, step = 1.000000, k = 112, ls = 1

Iteration 113:
  fx = 19644759.680134, xnorm = 141044.325394, gnorm = 2.923178, step = 1.000000, k = 113, ls = 1

Iteration 114:
  fx = 19644748.687705, xnorm = 141044.000170, gnorm = 2.878249, step = 1.000000, k = 114, ls = 1

Iteration 115:
  fx = 19644739.153086, xnorm = 141045.160075, gnorm = 2.830596, step = 1.000000, k = 115, ls = 1

Iteration 116:
  fx = 19644729.947724, xnorm = 141044.787943, gnorm = 2.468117, step = 1.000000, k = 116, ls = 1

Iteration 117:
  fx = 19644719.736514, xnorm = 141045.515515, gnorm = 2.185377, step = 1.000000, k = 117, ls = 1

Iteration 118:
  fx = 19644708.786184, xnorm = 141045.324937, gnorm = 2.203510, step = 1.000000, k = 118, ls = 1

Iteration 119:
  fx = 19644702.578301, xnorm = 141046.469059, gnorm = 2.604891, step = 1.000000, k = 119, ls = 1

Iteration 120:
  fx = 19644698.475381, xnorm = 141045.841024, gnorm = 2.724335, step = 1.000000, k = 120, ls = 1

Iteration 121:
  fx = 19644693.220169, xnorm = 141046.606666, gnorm = 2.402665, step = 1.000000, k = 121, ls = 1

Iteration 122:
  fx = 19644687.818987, xnorm = 141046.212669, gnorm = 1.993935, step = 1.000000, k = 122, ls = 1

Iteration 123:
  fx = 19644682.596253, xnorm = 141046.584436, gnorm = 1.678809, step = 1.000000, k = 123, ls = 1

Iteration 124:
  fx = 19644674.284326, xnorm = 141046.605784, gnorm = 1.561608, step = 1.000000, k = 124, ls = 1

Iteration 125:
  fx = 19644667.175717, xnorm = 141046.942135, gnorm = 1.535183, step = 1.000000, k = 125, ls = 1

Iteration 126:
  fx = 19644661.604728, xnorm = 141046.842505, gnorm = 1.484989, step = 1.000000, k = 126, ls = 1

Iteration 127:
  fx = 19644657.416867, xnorm = 141047.276729, gnorm = 1.604037, step = 1.000000, k = 127, ls = 1

Iteration 128:
  fx = 19644654.295803, xnorm = 141047.002335, gnorm = 1.603422, step = 1.000000, k = 128, ls = 1

Iteration 129:
  fx = 19644651.442649, xnorm = 141047.467161, gnorm = 1.455218, step = 1.000000, k = 129, ls = 1

Iteration 130:
  fx = 19644648.541117, xnorm = 141047.336586, gnorm = 1.261218, step = 1.000000, k = 130, ls = 1