% file:      	sparse_recovery.m, (c) Matthew Roughan, Tue Aug 22 2006
%
clear;
path('/home/mroughan/src/matlab/Michael_Saunders_Standford/', path);
path('/home/mroughan/src/matlab/NUMERICAL_ROUTINES/', path);

N = 3000;
x = (1:N)/N;
f = 3;
y = sin(2*pi*f*x);
y(floor(N/2.8)) = y(floor(N/2.8)) + 1;

figure(1)
plot(y,'b', 'linewidth', 3);
set(gca, 'linewidth', 3, 'xtick', [], 'ytick', []);
% axis off
set(gcf, 'PaperUnits', 'centimeters', 'PaperPosition', [0 0 25 10])
print('-depsc', sprintf('Plots/sparse_recovery.eps', i));

N = 30;
% N = 5;
x = (0:N-1)/N;
f_0 = 3;
y = sin(2*pi*f_0*x);
k_0 = floor(N/1.9);
y(k_0) = y(k_0) + 1;

figure(2)
hold off
plot([x; x], [zeros(size(y)); y],'b', 'linewidth', 3);
hold on
plot(x, y, 'bo', 'linewidth', 4);
set(gca, 'linewidth', 3, 'fontsize', 18);
set(gcf, 'PaperUnits', 'centimeters', 'PaperPosition', [0 0 28 10])
print('-depsc', sprintf('Plots/sparse_recovery_2.eps', i));


z = log10(abs(fftshift(fft(y))).^2);
figure(3)
plot(z);


%%% perform the sparse recovery via minimizing the l^1 norm
%%%%   alpha(1:N) are the coefficients of the deltas
%%%%   alpha(N+1:2*N) are the coefficeients of the sines
%%%%  
%%%%   equations
%%%%      y(k) = sum_{j=1}^N alpha(j) delta(k-j) + sum_{i=N+1}^{2*N}  alpha(j) exp(2*pi*sqrt{-1}*(i-N-1)*(j-1)/N)
%%%%           = alpha(k) + sum_{j=N+1}^{2*N}  alpha(j) exp(2*pi*i*(k-1)*(j-N-1)/N)
%%%%
%%%%   and
%%%%      -epsilon(i) <= alpha(i) <= epsilon(i)
%%%%
%%%%   we have extended our variables to be (alpha, epsilon)
%%%%   optimization criteria is 
%%%%          min sum_{i=1}^{2*N} epsilon(i)
%%%%
%%%%   4N variables
%%%%   6N constraints
%%%%
%%%% actually, cause can't do complex arithmetic in opt. do cos and sin
%%%%    6N variables
%%%%
f = [zeros(1,3*N) ones(1,3*N)];

[I,J] = meshgrid(0:N-1);
% W = exp(2*pi*sqrt(-1)*I.*J/N);
%   complex case not available
W1 = sin(2*pi*I.*J/N);
W2 = cos(2*pi*I.*J/N);
Aeq = [eye(N) W1 W2 zeros(N,3*N)];
beq = y';      

A1 = [eye(3*N) -eye(3*N)];
b1 = zeros(3*N,1);

A2 = [-eye(3*N) -eye(3*N)];
b2 = zeros(3*N,1);

A = [A1; A2];
b = [b1; b2];

[z, fval, exitflag, output, lamba] = linprog(f, A, b, Aeq, beq);

alpha = z(1:3*N);
deltas = alpha(1:N);
sines = alpha(N+1:2*N);
coses = alpha(2*N+1:3*N);

%%% but use the symmetry of the sine functions to remove redundant bits
k = ceil(N/2);
sines(2:end) = sines(2:end) - flipud(sines(2:end));

figure(5)
subplot(2,1,1)
hold off
plot([x; x], [zeros(size(deltas')); deltas'],'b', 'linewidth', 3);
hold on
plot(x, deltas', 'bo', 'linewidth', 4);
set(gca, 'linewidth', 3, 'fontsize', 18, 'ylim', [-1 1.6]);
ylabel('spike basis');

subplot(2,1,2)
hold off
plot([0:N/2-1; 0:N/2-1], [zeros(size(sines(1:N/2)')); sines(1:N/2)'],'b', 'linewidth', 3);
hold on
p1 = plot(0:N/2-1, sines(1:N/2)', 'bo', 'linewidth', 4);
p2 = plot(0:N/2-1, coses(1:N/2)', 'rx', 'linewidth', 2);
legend([p1 p2], 'sine', 'cosine');
set(gca, 'linewidth', 3, 'fontsize', 18, 'ylim', [-1 1.6]);
ylabel('Fourier basis');
set(gcf, 'PaperUnits', 'centimeters', 'PaperPosition', [0 0 28 16])
print('-depsc', sprintf('Plots/sparse_recovery_3.eps', i));



%% some debugging
actual_deltas = zeros(N,1);
actual_deltas(k_0) = 1;
actual_sines = zeros(N,1);
actual_sines(f_0+1) = 1;
actual_coses = zeros(N,1);
actual_alphas = [actual_deltas; actual_sines; actual_coses];

actual_solution = [actual_alphas; abs(actual_alphas)];
find(A*actual_solution > 0)
find(abs(Aeq*actual_solution-beq) > 0.001)
