function [mags,est,yh, sses, mses, ems, spars] = dwi_ebp( y,tb, est, xs, mags, xfinder ,nits, y0, est0, mags0)
% Elastic Basis Pursuit algorithm used for fitting DWI data (previously called generalized LH)
% This version computes validation errors
% Inputs:
%  - y, N x 1 vector: signal
%  - tb,, scalar: L1 norm of the true weights
%  - est, D x K_0 matrix: initial estimate of theta (the parameters)
%  - mags, 1 x K_0 vector: initial estimate of w ( the weights )
%  - xfinder, function handle: a function implementing the oracle tau, takes r (residual) as input
%  - nits, integer: number of iterations M
%  - y0, N x 1 vector: ground truth signal used for validation
%  - est0, D x K matrix: ground truth parameters
%  - mags0, 1 x K : ground truth weights
% Outputs:
%  - mags, 1 x K_1 vector: final estimate of weights w
%  - est, D x K_1 matrix, final parameter estimate theta
%  - yh, (N+1) x 1 vector: first N entries are predicted signal, last is L1 norm of w
%  - sses, 1 x M vector: sum of squares of residuals in each iteration
%  - mses, 1 x M vector: validated squared prediction error for each iteration
%  - ems, 1 x M vector:  validated earthmover distance for each iteration
%  - spars, 1 x M vector:  sparsity in each iteration
n=size(y,1);
r = y - xs(1:n,:) * mags;
[xnew,estnew,corr] = xfinder(r);
count=0;
ksize=size(mags,2);
oldyh = zeros(n+1,1);
sses = zeros(1,nits); mses = sses; ems = sses; spars=sses;
%%
for iter = 1:nits;
    %%
    est =[est,estnew];
    xs = [xs,[xnew;1]];
    mags = [mags;0];
    ksize=size(mags,1);
    ssp = lsqnonneg(xs, [y;tb]);
    xs = xs(:, ssp > 0);
    est = est(:, ssp > 0);
    mags = ssp(ssp > 0);
    yh = xs*mags;
    r=y-yh(1:n);
    spars(iter) = ksize;
    sses(iter) = norm([y;tb]-yh)^2 ;
    mses(iter) = norm(yh(1:n) - y0)^2;
    ems(iter) = symearthmover2(est0,mags0,est,mags,1e-5);
    [xnew,estnew,corr] = xfinder(r);
    oldyh=yh;
end
end

