function [ mats ] = PFR_mul_MakeMats( X, d, W, Query, varargin )  %
%PFR_MUL_MAKEMATS construct matrices for PFRank

ip = inputParser;
ip.FunctionName = 'PFR_mul_MakeMats';
ip.addRequired('X', @isnumeric);
ip.addRequired('d', @isnumeric);
ip.addRequired('W', @isnumeric);
ip.addRequired('Quary', @isnumeric);   %

ip.addOptional('T', [], @isnumeric);
ip.addOptional('use_parfor', false, @islogical);
ip.addOptional('use_sparse', false, @islogical);

ip.parse(X, d, W, Query, varargin{:});  %
par = ip.Results;

times = [];

%************************************************
% preparation
%************************************************
ts_prepare = cputime;

[n] = size(X,2);

% ensure W is not self-connected
W(sub2ind([n n], 1:n, 1:n)) = 0;

L  = diag(sum(W)) - W;

times.prepare = cputime - ts_prepare;

if ~isempty(par.T)
    Tcache = par.T;
    times.tangent_space = 0;
else
    ts_tangent_space = cputime;
    Tcache = TSEstimate(X', d, 'W', W);
    times.tangent_space = cputime - ts_tangent_space;
end

%************************************************
% matrix construction
%************************************************
ts_matrix = cputime;

if par.use_sparse
    mat_fun = @sparse;
else
    mat_fun = @zeros;
end

G     = mat_fun(d*n, d*n);
beta  = mat_fun(d*n, n);
B     = mat_fun(d*n, d*n);
D     = mat_fun(d*n, d*n);  %
H     = mat_fun(d*n, 1);    %
Id    = eye(d);

%dont understande
if par.use_parfor
    Gc     = cell(n,1);
    betac  = cell(n,1);
    Bc     = cell(n,1);
    
    parfor i = 1:n
        Gbk = zeros(d);
        bbk = zeros(d,n);
        Bbk = zeros(d);
        Xi  = X(:,i);
        Ti  = Tcache(:,:,i);
        
        nei = find(W(:,i))'; % sparse column access faster than row access
        Bcc = cell(size(nei));
        
        for ij = 1:length(nei)
            j = nei(ij);
            Mt1 = Ti'*(X(:,j)-Xi);
            Qij = Ti'*Tcache(:,:,j);
            Mt2 = zeros(d, n); Mt2(:,i) = -W(i,j)*Mt1; Mt2(:,j) = W(i,j)*Mt1;
            
            Gbk = Gbk + W(i,j)*(Mt1*Mt1');
            bbk = bbk + Mt2;
            Bbk = Bbk + W(i,j) * (Qij*Qij' + Id);
            Bcc{ij} = -2*W(i,j)*Qij;
        end
        
        Gc{i} = Gbk;
        betac{i} = bbk;
        Bc{i} = [Bbk Bcc];
    end
    
    for i = 1:n
        idx = (i-1)*d+1:i*d;
        nei = find(W(i,:));
        Bcc = Bc{i};
        for ij = 1:length(nei)
            j = nei(ij);
            jdx = (j-1)*d+1:j*d;
            B(idx,jdx) = Bcc{ij+1};
        end
        G(idx,idx)   = Gc{i};
        beta(idx,:)  = betac{i};
        B(idx,idx)   = Bcc{1};
    end
else
    for i = 1:n
        idx = (i-1)*d+1:i*d;
        Gbk = zeros(d);
        bbk = zeros(d,n);
        Bbk = zeros(d);
        Xi  = X(:,i);
        Ti  = Tcache(:,:,i);
        
        for j = find(W(i,:))
            jdx = (j-1)*d+1:j*d;
            Mt1 = Ti'*(X(:,j)-Xi);
            Qij = Ti'*Tcache(:,:,j);
            Mt2 = zeros(d, n); Mt2(:,i) = -W(i,j)*Mt1; Mt2(:,j) = W(i,j)*Mt1;
            
            Gbk = Gbk + W(i,j)*(Mt1*Mt1');
            bbk = bbk + Mt2;
            Bbk = Bbk + W(i,j) * (Qij*Qij' + Id);
            B(idx,jdx) = -2*W(i,j)*Qij;
        end
        
        G(idx,idx)   = Gbk;
        beta(idx,:)  = bbk;
        B(idx,idx)   = Bbk;
    end
end
%////////////////////////////////////////////////////////
for j = 1:n
    Dbk = zeros(d);
    Hbk = zeros(d, 1);
    Xj  = X(:,j);
    Tj  = Tcache(:,:,j);
    for t = find(Query')
        jdx = (j-1)*d+1:j*d;
        if (W(j,t)>0)
            Dbk = Id;
            Hbk = Hbk + Tj'*(X(:,t)-Xj)*Query(t);
            break;
        end;
    end;
    D(jdx, jdx) = Dbk;
    H(jdx) = Hbk;
end;        
%///////////////////////////////////////////////////////
times.matrix = cputime - ts_matrix;

mats = [];
mats.L = L;
mats.G = G;
mats.beta = beta;
mats.B = B;
mats.T = Tcache;
mats.W = W;
mats.D = D;      %
mats.H = H;      %
mats.times = times;
mats.d = d;
mats.n = n;

end
