% This is the GeLMA algorithm to solve an underdetermined system Ax =y
% x is sparse
% A is random
% The algorithm is used in 
% A differential equations approach to l_ 1-minimization with applications to array imaging
% by M.Moscoso, A.Novikov, G.Papanicolaou and L.Ryzhik 

n=200;% n columns, unknowns
delta=0.6;
spa=0.06;% sprsity parameter
m=round(delta*n); % m rows (measurements)
NN=linspace(1,n,n);

% Define matrix
A=randn (m,n)/sqrt(m);% we assume each coeff in [-1/sqrt(m),1/sqrt(m)]

%Define solution xx
k= round(sign (rand(n,1)-1+spa)/2+0.5);% sparsifier
xx=k.*(rand(n,1)-1/2);
koko=max(size(find(k>0)));
nok=1-k;
Nind=find(nok>0);

%Define right-hand side y
y=A*xx;
tau=mean(abs(A'*y));
dt=1/norm(A);

%Initialization
x=xx*0;
z=A*x;
ww=x;
LastRelativeError=0.1;
toler=0.000001;


while (LastRelativeError>toler)

znew=(y-A*x)*dt+z;
ww=A'*(y-A*x +z)*dt+x;
z=znew;
x=eta(ww,tau*dt);

figure(7)
  plot(NN,x,'k*',NN,xx,'g.')
  Title('Dots are true solution. Stars are approximation')
pause(0.1)

LastRelativeError=mean (abs(x-xx))/mean(abs(xx));

end


 
