%Message passing algorithm
%sf_g - first graph (Fly), selected proteins of ineterest
%sf_h - second graph (Yeast), selected proteins of ineterest
%sf_c - InParanoid assignment matrix (Fly-Yeast),sf_c(i,j)==1 if protein i
%(Fly) is a potential orthologus of prtoein j (Yeast)
%sf_inds_1 - index of selected proteins of interest in the entire list of
%proteins: Fly
%sf_inds_1 - index of selected proteins of interest in the entire list of
%proteins: Yeast
%sf_simple_match - all trivial assignments (all proteins without selected
%proteins)
function T=mp_algo(sf_g,sf_h,sf_c,sf_inds_1,sf_inds_2,sf_simple_match)
%%%%%%%%%%%%%%%%%%%DATA LOADING%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
G=sparse(load(sf_g));
H=sparse(load(sf_h));
C=sparse(load(sf_c));
G=G-diag(diag(G));
H=H-diag(diag(H));
Gt=G;Ht=H;
Gi=Gt;
Hi=Ht;
Ci=C;
cur_clust=0;clustG=zeros(size(Ci,1),1);clustF=zeros(size(Ci,2),1);clustSize=[];
for i=1:size(Ci,1)%length(ind_1)%size(Ci,1)
    if (clustG(i)==0)
        sel=find(Ci(i,:)>0);
        selg=find(full(sum(Ci(:,sel),2))>0);
        self=find(full(sum(Ci(selg,:),1))>0);
        cur_clust=cur_clust+1; 
        clustG(selg)=cur_clust;
        clustF(self)=cur_clust;
        clustSize(cur_clust)=length(selg)+length(self);
    end
end
%Matrix of Inparanoid cluster connectivity
Gcu=unique(clustG);Fcu=unique(clustF);
Inp_adj=zeros(length(Gcu),length(Fcu));
if 1
for i=1:length(Gcu)
    incG=clustG(find(sum(Gi(find(clustG==Gcu(i)),:),1)));
    for j=1:length(Fcu)
        %two clusters are connected if there is a coonection in G and in H
        Inp_adj(i,j)=max(max(Gi(find(clustG==Gcu(i)),find(clustG==Gcu(j)))))*max(max(Hi(find(clustF==Fcu(i)),find(clustF==Fcu(j)))))>0.25;
    end
    disp(i);
end
end

for i=1:length(Gcu)
    Inp_adj(i,i)=0;
end
%message passing algorithm
%binary variables
%association order is the following if there are more proteins of the first
%class then all possible combinations versus all possible permutations
%otherwise first class order is fixed and it's all combination*permutation
%of the second one
ind_1=1:length(clustG);ind_2=1:length(clustF);
Inp_adj_s=Inp_adj;
%kill all the stuff with unique possible assignement
ind_sel=1:size(Inp_adj,1);
Inp_adj(clustSize==2,:)=[];Inp_adj(:,clustSize==2)=[];
clustSizeF=clustSize(clustSize>2);
ind_sel=ind_sel(clustSize>2);
    for i=1:length(Gcu)
        clust_i1=unique(ind_1(clustG==Gcu(i)));
        clust_i2=unique(ind_2(clustF==Gcu(i)));
        l=1;
        for f=max(length(clust_i1),length(clust_i2))-min(length(clust_i1),length(clust_i2))+1:max(length(clust_i1),length(clust_i2))
            l=l*f;
        end
        V{i}=zeros(1,l);
        DefV(i)=0;
        if (l==1)
            V{i}=1;
            DefV(i)=1;
        end
    end
f=find(sum(Inp_adj)==0);
ind_sel_f=ind_sel(f);
backw=[];backw=[backw ind_sel(f)];BackInd=[];BackVert=[];
while ~isempty(Inp_adj)
    size(Inp_adj)
    for j=1:length(f)
        if (j==51)
               dbg=1;
        end
        finced=Gcu(find(Inp_adj_s(ind_sel_f(j),:).*DefV));
        %possible associations in cluster ind_sel_f(j)
        clust_i1_j=unique(ind_1(clustG==ind_sel_f(j)));
        clust_i2_j=unique(ind_2(clustF==ind_sel_f(j)));
        l_j=min(length(clust_i1_j),length(clust_i2_j));
        com_i1_j=combntns(clust_i1_j,l_j);
        com_i2_j=combntns_perm(clust_i2_j,l_j);
        ind_rep=1:size(com_i2_j,1);ind_rep=repmat(ind_rep,size(com_i1_j,1),1);ind_rep=ind_rep(:);
        com_i1_i2_j=[repmat(com_i1_j,size(com_i2_j,1),1) com_i2_j(ind_rep,:)]; 
        G_ind_j=repmat(com_i1_j,size(com_i2_j,1),1);G_ind_j=G_ind_j(:);
        H_ind_j=com_i2_j(ind_rep,:);H_ind_j=H_ind_j(:);
        %each accociation may produce intern conserved interactions as well
        CInt=0;
        %sum over size of association in G
        for h=1:l_j
            CIntt=G(G_ind_j((h-1)*size(com_i1_i2_j,1)+1:h*size(com_i1_i2_j,1)),G_ind_j).*H(H_ind_j((h-1)*size(com_i1_i2_j,1)+1:h*size(com_i1_i2_j,1)),H_ind_j)>0.26;
            CInt=CInt+CIntt;
        end
        %sum over size of association in H
        CIntp=CInt;CInt=0;
        for h=1:l_j
            CInt=CInt+CIntp(:,(h-1)*size(com_i1_i2_j,1)+1:h*size(com_i1_i2_j,1));
        end
        V{ind_sel_f(j)}=diag(CInt)';
        diag(CInt)';
        for k=1:length(finced)
            %possible associations in cluster finced(k)
            clust_i1_k=unique(ind_1(clustG==finced(k)));
            clust_i2_k=unique(ind_2(clustF==finced(k)));
            l_k=min(length(clust_i1_k),length(clust_i2_k));
            com_i1_k=combntns(clust_i1_k,l_k);
            com_i2_k=combntns_perm(clust_i2_k,l_k);
            ind_rep=1:size(com_i2_k,1);ind_rep=repmat(ind_rep,size(com_i1_k,1),1);ind_rep=ind_rep(:);
            com_i1_i2_k=[repmat(com_i1_k,size(com_i2_k,1),1) com_i2_k(ind_rep,:)]; 
            G_ind_k=repmat(com_i1_k,size(com_i2_k,1),1);G_ind_k=G_ind_k(:);
            H_ind_k=com_i2_k(ind_rep,:);H_ind_k=H_ind_k(:);
            %conserved interaction between pairs
            CInt=0;
            %sum over size of accociation in G
            for h=1:l_j
                CIntt=G(G_ind_j((h-1)*size(com_i1_i2_j,1)+1:h*size(com_i1_i2_j,1)),G_ind_k).*H(H_ind_j((h-1)*size(com_i1_i2_j,1)+1:h*size(com_i1_i2_j,1)),H_ind_k)>0.26;
                CInt=CInt+CIntt;
            end
            %sum over size of accociation in H
            CIntp=CInt;CInt=0;
            for h=1:l_k
                CInt=CInt+CIntp(:,(h-1)*size(com_i1_i2_k,1)+1:h*size(com_i1_i2_k,1));
            end
            %now CInt should have size size(com_i1_i2_j,1)*size(com_i1_i2_k,1)
    		v_add=CInt.*(repmat(V{finced(k)},length(V{ind_sel_f(j)}),1)+1); %+1 because we add conserved interactions on the new iteration
            %subraph in PPI alignment
            v_add=CInt+(repmat(V{finced(k)},length(V{ind_sel_f(j)}),1)); 
    		[v_add l_add]=max(v_add',[],1);
    		V{ind_sel_f(j)}=V{ind_sel_f(j)}+v_add;
    		%l_add(v_add==0)=-1;
    		BackInd{ind_sel_f(j)}(:,k)=l_add';
        end
        BackVert{ind_sel_f(j)}=finced;
        DefV(ind_sel_f(j))=1;
    end
    Inp_adj(f,:)=[];
    Inp_adj(:,f)=[];
    ind_sel(f)=[];
    f=find(sum(Inp_adj)<2);
    if ~isempty(Inp_adj)
        backw=[backw ind_sel(f)];
        ind_sel_f=ind_sel(f);
    end
end
SetV=zeros(size(DefV));%will be used for backward recursion
%backward recursion
for j=length(backw):-1:1
    finced=Gcu(find(Inp_adj_s(backw(j),:).*SetV));
    if isempty(finced)
        [Asv(backw(j)) As(backw(j))]=max(V{backw(j)});
    else
        if (length(finced)>1)
            disp('Initial graph is not a tree!!!!');
            return;
        end
        k=find(BackVert{finced}==backw(j));
        l=BackInd{finced}(As(finced),k);
        Asv(backw(j))=V{backw(j)}(l);
        As(backw(j))=l;
    end
    SetV(backw(j))=1;
end

%reconstruction of the optimal alignment
P_mp=[];InpClust=[];
for i=1:length(Gcu)
        clust_i1_j=unique(ind_1(clustG==Gcu(i)));
        clust_i2_j=unique(ind_2(clustF==Gcu(i)));
         if length(clust_i1_j)+length(clust_i2_j)<3
            P_mp=[P_mp; clust_i1_j clust_i2_j];
            continue;
         end
        l_j=min(length(clust_i1_j),length(clust_i2_j));
        com_i1_j=combntns(clust_i1_j,l_j);
        com_i2_j=combntns_perm(clust_i2_j,l_j);
        ind_rep=1:size(com_i2_j,1);ind_rep=repmat(ind_rep,size(com_i1_j,1),1);ind_rep=ind_rep(:);
        com_i1_i2_j=[repmat(com_i1_j,size(com_i2_j,1),1) com_i2_j(ind_rep,:)]; 
        opt_as=reshape(com_i1_i2_j(As(Gcu(i)),:),l_j,2);
        P_mp=[P_mp; opt_as];
end

ind_1_sel=load(sf_inds_1);
ind_2_sel=load(sf_inds_2);
T=load(sf_simple_match);
T=[T;[ind_1_sel(P_mp(:,1)) ind_2_sel(P_mp(:,2))]];

Ts=T(end,:);

for j=size(T,1):-1:1
   if ~ismember(T(j,1),Ts(:,1)) &&  ~ismember(T(j,2),Ts(:,2))
       Ts=[Ts;T(j,:)];
   end
end
T=Ts;
%all combinations where order is taking into account
function res=combntns_perm(x,m)
res=[];
if (m==1)
    if size(x,1)>size(x,2)
        res=x;
    else
        res=x';
    end
    return;
end
if length(x)<m
    if size(x,1)>size(x,2)
        res=x';
    else
        res=x;
    end
    return;
end

for i=1:length(x)
    xt=x;xt(i)=[];
    res_i=combntns_perm(xt,m-1);
    res_i=[ones(size(res_i,1),1)*x(i) res_i];
    res=[res;res_i];
end
