% Optical Imaging and Spectroscopy
%
% David J. Brady
% Duke University
% www.opticalimaging.org
%
% Figure 10.45
%
%
% This code modifies
% l1eq_example.m
%
% distributed at 
%
%
% Test out l1eq code (l1 minimization with equality constraints).
%
% Written by: Justin Romberg, Caltech
% Email: jrom@acm.caltech.edu
% Created: October 2005
%
% modified by DB July 2008 to compare global vs compact binary kernel
%
% put key subdirectories in path if not already there
close all;
figure(1);set(gcf,'color','white'); colormap gray;
%path(path, '.\Optimization');
%path(path, '.\Data');

% load random states for repeatable experiments
% load RandomStates
% rand('state', rand_state);
% randn('state', randn_state);

% signal length
N = 1024;
% number of spikes in the signal
T = 30;
% number of observations to make
K = 128;
% number of singular values to keep

% random +/- 1 signal
x = zeros(N,1);
q = randperm(N);
x(q(1:T)) = rand(T,1);

% measurement matrix
disp('Creating measurment matrix...');
A1=round(rand(K,N))/K;
[a1 b1 c1]=svd(A1);
A=a1(:,1:120)'*A1;
subplot(3,1,1);semilogy(diag(b1),'-k'); hold on;
title('(a) singular value spectra');
axis([-1 128 0.01 10]);
disp('Done.');

% observations
y = A*x;
figure(2);subplot(2,1,1);plot(y);hold on;
figure(1);
% initial guess = min energy
x0 =c1(:,1:120)*pinv(b1(1:120,1:120))*y;

% solve the LP
tic
xp = l1eq_pd(x0, A, [], y, 1e-3);
toc
subplot(3,1,2);plot([x'; (xp+1)']','-k');hold on;
title('(b) noise free reconstructions');
axis([-1 1024 -.5 3.2]);
%
% now try compact sampling kernel per pixel
%
% 8 measurement times, each with 16 measurements
% each pixel is assigned to 1 of the 16 detectors in each time step
%
A2=[(1/8)*ones(1,1024);zeros(15,1024)];
for pip=1:1024
    A2(:,pip)=A2(randperm(16),pip);
end
for pop=1:7
    Ap=[(1/8)*ones(1,1024);zeros(15,1024)];
    for pip=1:1024
        Ap(:,pip)=Ap(randperm(16),pip);
    end
    A2=[A2;Ap];
end
[a2 b2 c2]=svd(A2);
A=a2(:,1:120)'*A2;
subplot(3,1,1);semilogy(diag(b2),'-k');
% observations
y = A*x;
figure(2);subplot(2,1,2);plot(y);hold on;
figure(1);
% initial guess = min energy
x0 = c2(:,1:120)*pinv(b2(1:120,1:120))*y;

% solve the LP
tic
xp = l1eq_pd(x0, A, [], y, 1e-3);
toc
subplot(3,1,2);plot(xp+2,'-k');
%
%
% again with noise
%
%
% observations
y = a1(:,1:120)'*(A1*x+0.01*randn(128,1));
figure(2);subplot(2,1,1);plot(y);
figure(1);
% initial guess = min energy
x0 =c1(:,1:120)*pinv(b1(1:120,1:120))*y;
A=a1(:,1:120)'*A1;

% solve the LP
tic
xp = l1eq_pd(x0, A, [], y, 1e-3);
toc
subplot(3,1,3);plot(xp,'-k');hold on;
title('(c) reconstructions with noise variance \sigma^2=10^{-4}');
axis([-1 1024 -.5 2.2]);
% observations
y = a2(:,1:120)'*(A2*x+0.01*randn(128,1));
figure(2);subplot(2,1,2);plot(y);
figure(1);
% initial guess = min energy
x0 =c2(:,1:120)*pinv(b2(1:120,1:120))*y;
A=a2(:,1:120)'*A2;

% solve the LP
tic
xp = l1eq_pd(x0, A, [], y, 1e-3);
toc
subplot(3,1,3);plot(xp+1,'-k');