function [outT, outP0, outFMat] = PRW_Tri(inX, inYc, inYa, inL, inPara, inLearned)
% this code is to do the Preferenced random walk on Tri-relational graph
% input: 
% inX: n by d
% inYc: n by k
% inYa: n by l
% inL: the number of training instances. (we assume the data and labels are
%    re-arranged, that is, the first l instances are labeled and
%    the remaining are unlabeled.
% inPara:
%       inPara.k    : the k nearest neighbor graph
%       inPara.lambda:probability for (a,c) layers jumping
%       inPara.beta1: probability for (x,c) layers junmping
%       inPara.beta2: probability for (x,a) layers junmping
%       inPara.eta:   parameter for p0 init
%       inPara.gamma: probability for random walk restart
%       inPara.alpha:
% inLearned: the learned cell, the asymmetric more accuarate a-a, a-c, c-a,
%            c-c correlation
% output:
% outT: the big transition probability matrix
% outP0: the init matrix
% outFMat: the learned matrix after the label propogation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% check the size of the input arguments
if size(inX, 1) ~= size(inYa, 1)
	error('Incompatible input size for X and Ya');
else
    n = size(inX, 1);
    l = size(inYa, 2);
    fprintf('The number of instance is %d\n', n);
    fprintf('The number of annotation is %d\n', l);
end

if size(inX, 1) ~= size(inYc, 1)
    error('Incompatible input size for X and Yc');
else
    k = size(inYc, 2);
    fprintf('The number of class label is %d\n', k);
end
%% parameter settings
k = inPara.k;
lambda = inPara.lambda;
beta1 = inPara.beta1;
beta2 = inPara.beta2;
eta = inPara.eta;
gamma = inPara.gamma;
alpha = inPara.alpha;
%% build the affinity matrices    
tr_idx = 1: inL;
X_tr = inX(tr_idx, :);
Yc_tr = inYc(tr_idx,:);
Ya_tr = inYa(tr_idx,:);
% build the pairwise affinity matrix for the data graph, 
% using KNN graph and selftune method
X = inX';
X = X - repmat(mean(X, 2), 1, size(X, 2));
[wxx Wn] = selftuning(X', k);
if(nargin < 6)
    % build the pairwise affinity matrix for the annotation graph, using
    % consine kernel
    waa = getCosineKernel(inYa);
    % build the pairwise affinity matrix for the class label graph using Sb
    wcc = getYcCorr(Yc_tr, X_tr);
    % use cosine similarity as wca
    wac = zeros(size(Ya_tr, 2), size(Yc_tr, 2));
    for i = 1: size(Ya_tr,2)
        for j = 1: size(Yc_tr, 2)
            vec_a = Ya_tr(:,i);
            vec_c = Yc_tr(:,j);
            wac(i,j) = vec_a'*vec_c./(norm(vec_a)* norm(vec_c));
        end
    end
    wca = wac';
else
    waa = inLearned.f_matrix_aa;
    wcc = inLearned.f_matrix_cc;
    wac = inLearned.f_matrix_ac;
    wca = inLearned.f_matrix_ca;
end
% set the diagnal element of the affinity matrix to be 0;
wxx = wxx -diag(diag(wxx));
wcc = wcc -diag(diag(wcc));
waa = waa -diag(diag(waa));
% build the pairwise affinity matrix across the subgraph using the
% interationship via the label assignment
wxc = inYc;
wcx = inYc';
wxa = inYa;
wax = inYa';

%% Now, we build transition probability matrix
dxc = sum(wxc, 2);
dxc_zero_idx = (dxc == 0);
dxc(dxc_zero_idx) = eps;
txc = beta1 * diag(1 ./ dxc) * wxc;

dcx = sum(wcx, 2);
dcx_zero_idx = (dcx == 0);
dcx(dcx_zero_idx) = eps;
tcx = beta1 * diag(1 ./ dcx) * wcx;

dxa = sum(wxa, 2);
dxa_zero_idx = (dxa == 0);
dxa(dxa_zero_idx) = eps;
txa = beta2 * diag(1 ./ dxa) * wxa;

dax = sum(wax, 2);
dax_zero_idx = (dax == 0);
dax(dax_zero_idx) = eps;
tax = beta2 * diag(1 ./ dax) * wax;

dac = sum(wac,2);
dac_zero_idx = (dac == 0);
dac(dac_zero_idx) = eps;
tac = lambda * diag(1 ./dac) * wac;

dca = sum(wca,2);
dca_zero_idx = (dca == 0);
dca(dca_zero_idx) = eps;
tca = lambda * diag(1./dca) * wca;

txx = diag(1 ./ sum(wxx, 2)) * wxx;
tcc = diag(1 ./ sum(wcc, 2)) * wcc;
taa = diag(1 ./ (sum(waa, 2)+eps)) * waa;

tmp1 = (~dxc_zero_idx).*(~dxa_zero_idx);
tmp1 = (tmp1 >0);
tmp2 = (~dcx_zero_idx).*(~dca_zero_idx);
tmp2 = (tmp2 >0);
tmp3 = (~dax_zero_idx).*(~dac_zero_idx);
tmp3 = (tmp3 >0);

tmp4 = double(~dxc_zero_idx) - tmp1;
tmp4 = (tmp4 >0); 
tmp5 = double(~dxa_zero_idx) - tmp1;
tmp5 = (tmp5 >0);

tmp6 = double(~dcx_zero_idx) - tmp2;
tmp6 = (tmp6 >0); 
tmp7 = double(~dca_zero_idx) - tmp2;
tmp7 = (tmp7 >0);

tmp8 = double(~dax_zero_idx) - tmp3;
tmp8 = (tmp8 >0); 
tmp9 = double(~dac_zero_idx) - tmp3;
tmp9 = (tmp9 >0);

txx(tmp4, :) = (1- beta1) * txx(tmp4,:);
txx(tmp5, :) = (1- beta2) * txx(tmp5,:);
txx(tmp1, :) = (1- beta1 - beta2) * txx(tmp1,:);

tcc(tmp6, :) = (1- beta1) * tcc(tmp6,:);
tcc(tmp7, :) = (1- lambda) * tcc(tmp7,:);
tcc(tmp2, :) = (1- beta1 - lambda) * tcc(tmp2,:);

taa(tmp8, :) = (1- beta2) * taa(tmp8, :);
taa(tmp9, :) = (1- lambda) * taa(tmp9, :);
taa(tmp3, :) = (1- beta2 - lambda) * taa(tmp3,:);
% got the big trainsition probability matrix, the summation of each row should be 1
outT = [txx, txc, txa; tcx, tcc, tca; tax, tac, taa];

%% Construct p0 for init vector
Y = [inYc, inYa];
x0 = Y * diag(1 ./ sum(Y, 1));
y0 = eye(size(Y, 2));
outP0 = [eta * x0; (1 - eta) * y0]; % the summation of each column should be 1
%% label propogation
if(alpha == 0)
    outFMat = gamma * inv(eye(size(outT)) - (1 - gamma) * outT') * outP0;
else

tmpMat = [eye(size(wcc, 1)), zeros(size(wcc), size(waa));...
    zeros(size(a_corr,1), size(c_corr, 1)), a_corr];
eig_a_max = max(eig(waa));
alpha = min(alpha, 1 / eig_a_max - eps);
outFMat = gamma * inv(eye(size(outT)) - (1 - gamma) * outT') * ...
    outP0 * inv(eye(size(tmpMat)) - alpha * tmpMat);
end
end
