% SEND:
% - SP: string of the SP string that exist in all PP.
% - PP: cell array of PPs to analyze
% - Id: Consistency Matrix (output of SP-alignement functiony 'sp_pp_aline')

% DO:
% 0. estimate number of PPs that are inconsistent with one another.
% 1. for the consistent PPs: merge when possible the shorter PPs into the longer ones.
% 2. arrange final list of the fused PPs.
% 3. Get new Id matrix.
% 4. The folowing step is implemneted by "Main.m" function: go back to step 0 till all PPs are inconsistent

% GET:
% List of the new fused PPs 

function [new_cut_PP num_abs] = species_counting_byId_merge(SP, PP, Id)


% 0. lowest bound of number of species - number of inconsistent PPs
N = length(PP);
num_abs = N - length(union(Id(:,1), Id(:,2)));

% 1. Find the "blocks" of the Id matrix (all pair that match to a given PP)
%    - identify first column blocks
split = diff(Id(:,1));
%    - end-block index
ind2 = find(split > 0); dim = size(Id); dim = dim(1); ind2 = [ind2; dim];
%    - start-block index
ind1 = ind2+1; ind1 = [1; ind1]; ind1(end) = [];
%    - arrange by Block Size
R = [ind1 ind2 ind2-ind1]; R = sortrows(R,3); ind1 = R(:,1); ind2 = R(:,2);

% 2. Algo: 
%    A. go over each block.
%    B. for each single PP in the block, get all PPs that match to it [in Id(block,2)].
%    C. tag the PPs to merge and any PP that was left un-merged.
%       These constitute the new list of PPs to analyze in next iteration.
%    D. The number of sequences at the end that are left with no-match is the number of species


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%          STEP 1         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%   Identify the short PPs that can be absorbed/merged to longer ones     %%%%%%%%%%%% 

% go over each block
Nb = length(ind1);
% new PP list and the ones to avoid (after matching)
new_PP = [];
avoid = [];
tomerge = [];
k = 0;
for i = 1:Nb
        
    % get the block to treat => unify all of its PPs
    block = Id(ind1(i):ind2(i),:);
    union_block = union(block(:,1), block(:,2));
    Mb = length(union_block);
    
    % Skip PPs that has been already tagged as absorbed/merged with others (no point to analyze them again)
    skip = find( avoid == block(1,1) );
    if ~isempty(skip)
        continue;
    end
    
    % arrange PPs by their length
    for j = 1:Mb
        L(j) = length(PP{union_block(j)});
    end
    if length(L) <=2
        tmp = [union_block' L'];
    else
        tmp = [union_block L'];
    end
	tmp = sortrows(tmp,2);
    union_block = tmp(:,1);
    clear tmp j L
    
    % create consistency matrix of all possible matches of PPs belonging to the current block 
    % (PP index are in "union_block" => arranged by length)
    BlockMat = zeros(Mb,Mb);
    for j = 1:Mb
        % index - of all the other match of the PPs in the block
        x = Id(find(Id(:,1) == union_block(j)),:); x1 = union(x(:,1),x(:,2));
        y = Id(find(Id(:,2) == union_block(j)),:); y1 = union(y(:,1),y(:,2));
        z = union(x1,y1);
        
        [num id] = intersect(union_block,z);
        BlockMat(j,id) = 1;
    end      
    
    % go line by line (from the shortest to the longest).
    % - Compare each (short) PP with the longest it matches to and: 
    %   a. discard it if contained in the longer.
    %   b. otherwise, tag it to be merged with the longer. 
    % - Tag all the PP that should be kept and those that can be avoided
    count = 0; 
    keep = 0;
    for j = 1:Mb
        
        % get the longest possible to match to the current PP
        tmp = find(BlockMat(j,:) == 1); tmp = tmp(end);
        
        % check if the current PP is completely included in the longest one
        isincluded = findstr(PP{union_block(j)}, PP{union_block(tmp)});
        if isempty(isincluded)
            %display('merge by SP-alignement is needed');
            k = k + 1;
            tomerge(k,:) = [ union_block(j) union_block(tmp) ];
        end
        
        % check if merged to an old or new PP
        isalready = find (tmp == keep);
        if isempty(isalready)
            count = count + 1;
            keep(count) = tmp;
        end
        
    end
    
    % create new list of fused PPs
    tmp_new_PP = union_block(keep);
    new_PP = union(new_PP, tmp_new_PP);
    
    % create list of PPs have been already analyzed and can therefore be skipped in next "blocks" anlysis
    tmp_avoid = setdiff(union_block, tmp_new_PP);
    avoid = union(avoid, tmp_avoid);
    
    clear tmp* union* b* B* x* y* z*
end
clear ind*

% Get the new_PP sequences
new_cut_PP = PP(new_PP);


%%%%%%%%%%%%%%%%%%%%           STEP 2           %%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%   analyse the PP that require merging   %%%%%%%%%%%%%%%%%%%
count = 0;
keep = 0;
if ~isempty(tomerge)
    % Arrange the "tomerge" array by the second column (the longer protein)
    tomerge = sortrows(tomerge,2);

    % get the blocks of "tomerge"
    split = diff(tomerge(:,2));
    % end-block index
    ind2 = find(split > 0); dim = size(tomerge); dim = dim(1);
    ind2 = [ind2; dim];
    % start-block index
    ind1 = ind2+1;
    ind1 = [1; ind1];
    ind1(end) = [];

    % Go over the "tomerge" blocks array
    N = length(ind1);
    for i = 1:N
        
        block = tomerge(ind1(i):ind2(i), :);

        % 1) identify where is the longer PP in "new_PP"
        idx = block(1,2);
        pp_new = find (new_PP == idx);

        % 2) arrange the PP by length, starting from the short cndidate 
        %    - get the longest possible merge of these PP. 
        %    - Keep those that can not be merged
        Mb = size(block); Mb = Mb(1);
        for j = 1 : Mb
            L(j) = length(PP{block(j,1)});
        end
        % arrange by length
        tmp = [block(:,1) L'];
        tmp = sortrows(tmp,2);
        pp_merge = tmp(:,1);

        % merge process
        if Mb > 1
            pp_tmp = PP{block(1,1)};
            PPs{1} = pp_tmp;
            for j = 1:Mb-1
                % check consistency of the current PP, 'PP{block(j+1,1)}' with the updated merged PP 'pp_tmp'
                PPs{2} = PP{block(j+1,1)};
                tmp_Id = sp_pp_aline(SP, PPs);

                % if consistent - merge
                if ~isempty(tmp_Id)
                    pp_tmp = merge(pp_tmp, PP{block(j+1,1)}, SP);
                % if inconsistent 
                % - do not merge and keep seperately from the merged, if not already tagged as kept/new-pp
                else
                    %display('Can not merge -  - keep for next iteration');
                    iskept = find(keep ==  block(j+1,1));
                    isnewpp = find(new_PP ==  block(j+1,1));
                    if isempty(iskept) && isempty(isnewpp)
                        count = count +1;
                        keep(count) = block(j+1,1);
                    end
                end
            end
        else
            pp_tmp = PP{pp_merge};
        end

        % merge it with its "match"
        pp_tmp = merge(pp_tmp,PP{idx}, SP);

        % 3) put the merged PP in the right place in new PP
        new_cut_PP{pp_new} = pp_tmp;
        
        clear L Mb tmp* idx block pp_*
    end
    
    % 4) Keep those that could not have been merged.
    if count > 0
        add_PP = PP(keep);
        new_cut_PP = [new_cut_PP add_PP];
    end
end
