function [x,est,corr] = xfinder_random2( grid,rad,nits,newits, r )
% [x,est,corr] = xfinder_random2( grid,rad,nits,newits, r )
%  Oracle for EBP applied to DWI data
%  This function is not used directly, but rathe used to create function handles
%  e.g.
%  xfinder = @(r) xfinder_random2(grid,2,100,5,r);
% 
% Inputs:
%  - grid, 3xn matrix of measurement directions x_1,...,x_n
%  - rad, scalar: this should be the square root of the maximum possible diffusivity
%  - nits: number of random starting points
%  - newits: number of Newton iterations
%  - r, nx1 residual
%
% Outputs:
%  - x, nx1 kernel signal: the kernel signal with the highest inner product with r
%  - est, 3x1 vector: the diffusion-weighted fiber direction generating x
%  - corr =  x'*r

% generate a set of candidate fibers randomly
randcenters = rad.*2.*(rand(3,100.*nits)-1);
nms = sqrt(sum(randcenters.^2));
randcenters = randcenters(:,nms < rad);
randcenters = randcenters(:,1:nits);
xs = ste_tan0(randcenters,grid);
rvs = xs'*r;
corr = max(rvs);
x = xs(:,find(rvs==corr,1));
est = randcenters(:,find(rvs==corr,1));

candidates = randcenters;

% some functions which will be used for computing Newton directions
detsf = @(rd2s) rd2s(:,1).*(rd2s(:,4).*rd2s(:,6)-rd2s(:,5).^2) - rd2s(:,2).*(rd2s(:,2).*rd2s(:,6) - rd2s(:,5).*rd2s(:,3))...
        +rd2s(:,3).*(rd2s(:,5).*rd2s(:,2) - rd2s(:,4).*rd2s(:,3));
adjgsf = @(rd2s) [rd2s(:,4).*rd2s(:,6)-rd2s(:,5).^2,...
        rd2s(:,3).*rd2s(:,5)-rd2s(:,2).*rd2s(:,6),...
        rd2s(:,2).*rd2s(:,5)-rd2s(:,3).*rd2s(:,4),...
        rd2s(:,1).*rd2s(:,6)-rd2s(:,3).^2,...
        rd2s(:,2).*rd2s(:,3)-rd2s(:,1).*rd2s(:,5),...
        rd2s(:,1).*rd2s(:,4)-rd2s(:,2).^2];
adjsf = @(adjgs,dets,rd1s) (-1/2).*[adjgs(:,1).*rd1s(:,1)+adjgs(:,2).*rd1s(:,2)+adjgs(:,3).*rd1s(:,3),...
        adjgs(:,2).*rd1s(:,1)+adjgs(:,4).*rd1s(:,2)+adjgs(:,5).*rd1s(:,3),...
        adjgs(:,3).*rd1s(:,1)+adjgs(:,5).*rd1s(:,2)+adjgs(:,6).*rd1s(:,3)]./repmat(dets,1,3);
rd1sf = @(d1s) [d1s{1}'*r,d1s{2}'*r,d1s{3}'*r];
rd2sf = @(d2s) 1/2.*[d2s{1,1}'*r,d2s{1,2}'*r,d2s{1,3}'*r,d2s{2,2}'*r,d2s{2,3}'*r,d2s{3,3}'*r];
vstats = @(adjs,rd1s,rd2s) sum(adjs.*rd1s,2) + rd2s(:,1).*adjs(:,1).^2 ...
    +rd2s(:,4).*adjs(:,2).^2 +rd2s(:,6).*adjs(:,3).^2 +...
    2.*rd2s(:,2).*adjs(:,1).*adjs(:,2)+...
    2.*rd2s(:,3).*adjs(:,1).*adjs(:,3)+...
    2.*rd2s(:,5).*adjs(:,2).*adjs(:,3);


if size(candidates,2) ==0;
    x = ste_tan0(est,grid);
    return;
end

% compute kernels and keep the kernel (and fiber) with highest inner product
[ys,d1s,d2s] = ste_tan0(candidates,grid);
max(rvs);
best_rvs = ys'*r;
best_candidates = candidates;
%%
for iii=1:newits;
    %% compute Newton directions for the candidate fibers
    rd1s = rd1sf(d1s);
    rd2s = rd2sf(d2s);
    
    dets = detsf(rd2s);
    adjgs = adjgsf(rd2s);

    adjs = adjsf(adjgs,dets,rd1s);

    % move candidate fibers in Newton direction
    candidates = candidates + adjs';
    nms = sqrt(sum(candidates.^2));
    tempfilt = (nms > rad);
    if sum(tempfilt > 0);
        % rescale fibers if they exceed the maximum diffusivity  
        candidates(:,tempfilt) = candidates(:,tempfilt).*rad./repmat(nms(tempfilt),3,1);
    end
    [ys,d1s,d2s] = ste_tan0(candidates,grid);
    rvs = ys'*r;
    max(rvs);
    replfilt = (rvs > best_rvs);
    best_rvs(replfilt) = rvs(replfilt);
    best_candidates(:,replfilt) = candidates(:,replfilt);

    % update the best fiber found so far    
    est = best_candidates(:,find(best_rvs==max(best_rvs),1));
    corr=max(best_rvs);
end

x = ste_tan0(est,grid);

end

