% deconvolution_data.m
%
% Data matrix providing A, B, and Y.

n = 2048;  % Length of vector x
m = 256;  % Number of measurements
d = 16;  % Image dimension

% Knocked out sensor entries (these correspond to 0 in the Y vector).

knockouts = [1, 8, 17, 26, 43, 47,49, 51, 52, 54, 56, 61, 67, 83, 84,95,96, 104, 105, 106, 119, 121, 126, 136, 148, 155, 164, 171, 173, 175, 176, 184, 192, 193, 198, 202, 204, 207, 208, 224, 230, 238, 239, 240, 242, 244, 246, 249, 250, 256];

% Stacked version of the observed Y matrix. All knockout entries are 0.

Y = [0.0, 1.0, 0.986736, 0.950834, 0.890618, 0.802241, 0.730228, 0.0, 0.617109, 0.572867, 0.537012, 0.505473, 0.479585, 0.460734, 0.434283, 0.407227, 0.0, 0.943112, 0.93028, 0.900096, 0.850174, 0.776685, 0.715532, 0.661522, 0.616168, 0.0, 0.543681, 0.514492, 0.490327, 0.472822, 0.446446, 0.418632, 0.917611, 0.924769, 0.908196, 0.876034, 0.825046, 0.750903, 0.695771, 0.654471, 0.618582, 0.585544, 0.0, 0.532771, 0.511308, 0.495924, 0.0, 0.440252, 0.0, 0.942465, 0.0, 0.0, 0.822465, 0.0, 0.679012, 0.0, 0.608087, 0.58723, 0.569284, 0.550538, 0.0, 0.522534, 0.496656, 0.465714, 0.925014, 0.915113, 0.0, 0.848458, 0.79261, 0.703577, 0.641276, 0.599525, 0.571709, 0.555844, 0.547348, 0.543127, 0.537572, 0.534455, 0.511607, 0.479733, 0.901381, 0.889765, 0.0, 0.0, 0.762819, 0.671031, 0.607337, 0.565367, 0.539084, 0.526613, 0.524472, 0.529723, 0.538817, 0.54695, 0.0, 0.0, 0.870048, 0.86508, 0.836626, 0.794755, 0.731969, 0.639567, 0.575855, 0.0, 0.0, 0.0, 0.502627, 0.514456, 0.533332, 0.556715, 0.544042, 0.510148, 0.88082, 0.881872, 0.852725, 0.806534, 0.734776, 0.63639, 0.0, 0.524304, 0.0, 0.488821, 0.493562, 0.508894, 0.533496, 0.0, 0.5615, 0.526518, 0.861447, 0.863572, 0.834812, 0.790346, 0.719824, 0.620365, 0.551936, 0.0, 0.482313, 0.474154, 0.481509, 0.49987, 0.528422, 0.566848, 0.573385, 0.537663, 0.867501, 0.869613, 0.83766, 0.0, 0.71282, 0.611836, 0.54253, 0.497902, 0.472972, 0.46596, 0.0, 0.496203, 0.527174, 0.56849, 0.579279, 0.54319, 0.842554, 0.839698, 0.80547, 0.0, 0.686925, 0.589376, 0.522965, 0.480936, 0.458674, 0.454589, 0.0, 0.491243, 0.0, 0.566392, 0.0, 0.0, 0.859132, 0.850724, 0.811988, 0.766404, 0.69285, 0.593045, 0.525301, 0.0, 0.460709, 0.45775, 0.473277, 0.49783, 0.53043, 0.570177, 0.56848, 0.0, 0.0, 0.860428, 0.819414, 0.77263, 0.693422, 0.0, 0.525367, 0.483245, 0.462323, 0.0, 0.479338, 0.0, 0.534685, 0.567636, 0.0, 0.0, 0.863925, 0.858327, 0.818639, 0.769776, 0.690254, 0.59092, 0.524347, 0.483664, 0.464622, 0.465863, 0.487417, 0.512339, 0.537492, 0.557952, 0.544007, 0.0, 0.884122, 0.87514, 0.830235, 0.783669, 0.703779, 0.0, 0.534685, 0.493685, 0.47514, 0.477726, 0.50167, 0.523405, 0.537414, 0.0, 0.0, 0.0, 0.898131, 0.0, 0.831866, 0.0, 0.71376, 0.0, 0.542731, 0.501658, 0.0, 0.0, 0.513124, 0.526961, 0.535028, 0.542364, 0.523565, 0.0]';

% Now, let's construct the filter matrix.
%
% The translation is as follows, assuming a zero-indexed image (we
% must add one to the coordinates in Matlab as it does 1-indexing):
% from coordinate (i, j) in {0, ..., d-1}^2 in the image, we get to
% vector position
%
%  k = d * i + j.
%
% Thus from vector position k, image coordinates are
%
%   j = mod(k, d).
%   i = (k - j) / d.
%
% We construct row l, corresponding to position (i, j) in the image,
% of A so that for entry k corresponding to index pair (i', j') in the
% image, we have
%
%  A(l, k) = exp(-abs(i - i') / sigma - abs(j - j') / sigma)
%
% where sigma is the bandwidth of the filter / convolution.  Note that
% Matlab does 1-indexing, so we have to add 1 to l and k when we write
% A(l, k).

bandwidth = 5.0;

A = zeros(m, m);

for ll = 0:(m-1)
  j1 = mod(ll, d);
  i1 = floor((ll - j1) / d);
  for i2 = 0:(d-1)
    for j2 = 0:(d-1)
      kk = d * i2 + j2;
      A(ll+1, kk+1) = exp(-abs(i1 - i2) / bandwidth - abs(j1 - j2) / bandwidth);
    end
  end
  A(ll+1, :) = A(ll+1, :) / sum(A(ll+1, :));
end

A(A <= 1e-5) = 0;

A(knockouts, :) = 0;

% Now load the B matrix.
B = importdata("B_matrix_deconv_data.csv");

