%**************************************************************************
% WITH THE SIMPLE NETWORK EVOLUTIONARY MODEL, NO CONFIDENCE VALUES       *
% add the "true" leaves to the network formed by phylogeny, there will   *
% be different substitution probabilities from true leaves to noisy ones *
% In the framework, use the "true" leaves to replace the noisy ones, and *
% use the L() and C() of the "true" leaves to replace those of the noisy *
% ones (these new L and C should be calculated from the old ones). In the*
% traceback process,
% INPUT: tree -- the phylogenetic tree, represented as a matrix of size
%               #(nodes in the tree)x5. The nodes in the tree are numbered.
%               Each row in the matrix corresponds to a tree node:
%               [(idx of current node) (idx of its parent) (idx of its left
%               child) (idx of right child) (edge length to its parent)]
% 		The index of a parent node should be larger than both indices
%               of its children nodes.
%        nleaves - #(leaves in the tree)
%        leaf_net_noisy -- the input noisy networks for the leaf organisms.
%               It's a cell variable of size 1xnleaves, each cell entry has
%               the adjacency matrix of the corresponding leaf network.
%	 varargin - optional input. One can give entry-specific confidence 
%               values for the input noisy leaf networks in addition to the
%               DiffRate parameter. If this input is given, the DiffRate 
%               parameter will be modified to take into account these values.
%**************************************************************************

function [loglh allnets] = ProPhyC(tree, nleaves, leaf_net_noisy, varargin)

global Pii Prate DiffRate;
% main parameters:
% Pii - initial distribution of 0s and 1s
% Prate - the evolutionary rates (substitution matrix)
% DiffRate - the noise rate between the true leaf networks and the input noisy leaf networks
% Example of parameter setting:
% Pii = [0.9 0.1];
% Prate = [0.94 0.06; 0.05 0.95];
% DiffRate = [0.9 0.1; 0.3 0.7];

ntreenodes = size(tree,1);
ngenes = size(leaf_net_noisy{1},1);
dr_seq = cell(1,nleaves);
confi_seq = cell(1,nleaves);


% modify the DiffRate parameter if interaction-specific confidence values are given
k_inc10 = 1/8; % proportion to increase for P10
k_dec10 = 1/5; % proportion to decrease for P10
k_inc01 = 1/7; % proportion to increase for P01
k_dec01 = 1/8; % proportion to decrease for P01
if ~isempty(varargin) % refine with confidence values.
    net_confi = varargin{1};
    for l = 1:nleaves
        confi_seq{l} = matrix2seq(net_confi{l}); % take the absolute value!
        for s = 1:ngenes*ngenes
            dr_seq{l}{s} = DiffRate;
            p10 = DiffRate(2,1);
            p01 = DiffRate(1,2);
            if confi_seq{l}(s) >= 0.6 && confi_seq{l}(s) <= 0.9
                % confidence value (cv) for an edge, change p10
                % assume the values are between 0.6 and 0.9.
                % the higher the cv, the smaller p10 is.
                % when cv=0.6, p10=p10+p10*k_inc10
                % when cv=0.9, p10=p10-p10*k_dec10
                dr_seq{l}{s}(2,1) = p10 + p10*(1/3)*...
                    (-10*(k_inc10+k_dec10)*confi_seq{l}(s)+(9*k_inc10+6*k_dec10));
                
            elseif confi_seq{l}(s) > 0.9
                dr_seq{l}{s}(2,1) = p10 - p10*k_dec10;
                
            elseif confi_seq{l}(s) < 0.6 && confi_seq{l}(s) > 0
                dr_seq{l}{s}(2,1) = p10 + p10*k_inc10;
                
            elseif confi_seq{l}(s) >= -0.9 && confi_seq{l}(s) <= -0.6
                % change p01
                dr_seq{l}{s}(1,2) = p01 + p01*(1/3)*...
                    (10*(k_inc01+k_dec01)*confi_seq{l}(s)+(9*k_inc01+6*k_dec01));
                
            elseif confi_seq{l}(s) > -0.6 && confi_seq{l}(s) < 0
                dr_seq{l}{s}(1,2) = p01 + p01*k_inc01;
            elseif confi_seq{l}(s) < -0.9
                dr_seq{l}{s}(1,2) = p01 - p01*k_dec10;
            end
            dr_seq{l}{s}(2,2) = 1-dr_seq{l}{s}(2,1);
            dr_seq{l}{s}(1,1) = 1-dr_seq{l}{s}(1,2);
        end
    end
else % every entry in every leaf have the same diff_rate
    for l = 1:nleaves
        for s = 1:ngenes*ngenes
            dr_seq{l}{s} = DiffRate;
        end
    end
end


% take into account the edge lengths of the tree
max_length = max(tree(:,5));
max_length = int8(max_length);
P = cell(1,max_length+1);
P{1} = [1 0; 0 1]; % P{1} is actually P{0}
P{2} = Prate;
for i = 3:int8(max_length+1)
    P{i} = substi_prob_L(Prate, double(i-1));
end
edge_length = tree(:,5);
edge_length = edge_length+ones(ntreenodes,1);
edge_length = int8(edge_length);


allnets = cell(ntreenodes,1);
for i = 1:ntreenodes
    allnets{i} = zeros(ngenes);
end
loglh = 0;
for i = 1 : ngenes 
    for ii = 1:ngenes
        C = zeros(ntreenodes,2); 
        L = zeros(ntreenodes,2);
        C0 = zeros(nleaves,2); 
        L0 = zeros(nleaves,2);
        for j = 1:ntreenodes 
            if j >= 1 && j <= nleaves % for the leaves
                C0(j,1)=leaf_net_noisy{j}(i,ii)+1;
                C0(j,2)=leaf_net_noisy{j}(i,ii)+1;
                L0(j,1)=dr_seq{j}{i}(1,(leaf_net_noisy{j}(i,ii)+1));
                L0(j,2)=dr_seq{j}{i}(2,(leaf_net_noisy{j}(i,ii)+1));
                
                [L(j,1) C(j,1)] = max([P{edge_length(j)}(1,1)*L0(j,1)...
                    P{edge_length(j)}(1,2)*L0(j,2)]);
                [L(j,2) C(j,2)] = max([P{edge_length(j)}(2,1)*L0(j,1)...
                    P{edge_length(j)}(2,2)*L0(j,2)]);
                
            elseif j >= nleaves+1 && j <= ntreenodes-1 % for nonroot internal nodes:
                [L(j,1), C(j,1)] = max([P{edge_length(j)}(1,1)*L(tree(j,3),1)*L(tree(j,4),1)...
                    P{edge_length(j)}(1,2)*L(tree(j,3),2)*L(tree(j,4),2)]);
                [L(j,2), C(j,2)] = max([P{edge_length(j)}(2,1)*L(tree(j,3),1)*L(tree(j,4),1)...
                    P{edge_length(j)}(2,2)*L(tree(j,3),2)*L(tree(j,4),2)]);
            elseif j == ntreenodes
                [L(j,1), C(j,1)] = max([Pii(1)*L(tree(j,3),1)*L(tree(j,4),1) ...
                    Pii(2)*L(tree(j,3),2)*L(tree(j,4),2)]);
                % leave the second element as 0
            end
        end
        % traceback: from end to the beginning.
        allnets{ntreenodes}(i,ii)=C(ntreenodes,1);
        for j = ntreenodes-1:-1:1
            allnets{j}(i,ii) = C(j,allnets{tree(j,2)}(i,ii));
        end
        loglh = loglh+log2(L(ntreenodes,1));
    end
end


clear leaf_net_noisy;
for i = 1:ntreenodes
    allnets{i} = allnets{i}-1;
end



