%% new file for NIPS rebuttal

addpath('GridSphere/GridSphere');
load('realdata0.mat');
grid = bvecss{1,1};
rad  = 2;
centers  = sphlattice(10,rad);
xs       = ste_tan0(centers,grid);
biggrid  = biggridsphere(1000);
thetas   = (1:100)./100 .* 2*pi;
circgrid = [sin(thetas);cos(thetas);0.*thetas];
pp=size(centers,2);
% utility functions
applyToGivenRow = @(func, matrix) @(row) func(matrix(row, :));
newApplyToRows = @(func, matrix) arrayfun(applyToGivenRow(func, matrix), 1:size(matrix,1), 'UniformOutput', false)';
takeAll = @(x) reshape([x{:}], size(x{1},2), size(x,1))';
genericApplyToRows = @(func, matrix) takeAll(newApplyToRows(func, matrix)); % from http://stackoverflow.com/users/376454

rng(0);
% The number of fibers
true_k= 3;

% Partial volume fractions for the fibers
true_amp  = [.4;.3;.3];
tb        = sum(true_amp);
true_tamp = sum(true_amp);
% We simulate three fibers with slightly different directions
%In this parameterization we combine directionality of the fiber with
% its axial diffusivity
%
% If v is the parameter of the fiber, sqrt(norm(v))=b(lambda1-lambda2)
true_pos=[0.4657    1.2493    0.4667;
          -1.6289    0.6602   -1.1222;
          -0.8092    0.9064    1.2146];
%true_pos = [[1.52;-0.153;0.1],[.8;1.2;0.111],[.8;0.23;1.3]];
%true_pos = [ 0.4657    1.2493    0.4667];

%s Noise.
sigma=0.005;
noise=sigma;

%%
res_0 = zeros(1,100);
res_e = zeros(1,100);

%%
parfor ii=1:100;
    %%
    [data_orig,data_noiseless] = generate_data(true_pos,true_amp,grid,noise);
    y = data_orig; % signal with noise
    y0 = data_noiseless; % noiseless signal
    %% Fit NNLS model
    l1p=0.01;
    betaa_0 = lsqnonneg([xs;ones(1,pp)],[data_orig;1]);
    est_0 = centers(:,betaa_0 > 0);
    mags_0 = betaa_0(betaa_0 > 0);
    spar_0 = sum(betaa_0 > 0);
    xs_0 = xs(:,betaa_0 > 0);
    em_0 = symearthmover2(true_pos,true_amp,est_0,mags_0,1e-3);
    yh_0 = xs * betaa_0;
    mse_0 = norm(y-yh_0)^2;
    res_0(ii) = em_0;

    %% Fit EBP
    xfinder = @(r) xfinder_random2(grid,2,100,3,r);
    [mags_e, est_e, yh_e, sses_e, mses_e, ems_e, spars_e] = dwi_ebp(y, 1, zeros(3,1), ones(151,1), mean(y), xfinder, 1000, y0, true_pos, true_amp);
    spar_e = size(est_e,2);
    em_e = symearthmover2(true_pos,true_amp,est_e,mags_e,1e-3);
    res_e(ii) = em_e;
end

