% Syntax: [levels, n_fb, out]=hierarchy_levels(s1,nlevels,itop,ibot,beta, niter,weigth_len,levels);
% Sorts nodes in the network in NLEVELS levels to minimize the number of
% edges that go against the hierarchy, i.e. connecting a node from a lower level to this on 
% a higher or the same level.
% s1 - the adjacency matrix of a network
% nlevels - the number of hierarchical levels 
% OPTIONAL: (give [] for no top nodes) itop - indices of nodes to be placed at top level (level number nlevels)
% OPTIONAL: (give [] for no bottom nodes) ibot - indices of nodes to be placed at bottom level (level number 1)
% OPTIONAL: beta - inverse temperature of the Metropolis algorithm(DEFAULT: 10)
% OPTIONAL: niter - number of iterations (DEFAULT: 10000)
% OPTIONAL: levels - starting level array, NB. if levels is specified, itop
% and ibot cannot be used
% OPTIONAL: weight_len, punish long feedforward edges. (DEFAULT=0);
% NB, auto reg is removed at the beginning..
% K-K Yan adopted from Sergei Maslov. Ref: Ispolatov and Maslov, BMC
% Bioinformatics 2008

function [levels, n_fb,out]=hierarchy_levels(s1,nlevels,itop,ibot,beta,niter,weight_len,levels)

s1=sign(abs(s1-diag(diag(s1))));
%n_edges=full(sum(sum(s1)));
len1=length(s1);
ap=sign(sum(s1+s1'));%look at nodes participated in the network...
a_changeable=ap;

if (nargin < 5) beta=10; end;

if (nargin < 6) niter=10000; end;

if (nargin<7) weight_len=0; end;

if (nargin < 8)
        levels=zeros(len1,1);
        if isempty(itop)==0
                levels(itop)=nlevels;
                a_changeable(itop)=0;
        end;
        
        if isempty(ibot)==0
                levels(ibot)=1;
                a_changeable(ibot)=0;
        end;
        
        ip=find(a_changeable);
        levels(ip)=floor((nlevels).*rand(size(ip)))+1;

end        
        
        
i_ch=find(a_changeable)';
n_ch=length(i_ch);

[a,b]=find(s1);

n_fb_orig=sum(levels(a)<=levels(b));
tot_len_orig=sum((levels(a)>levels(b)).*(levels(a)-levels(b)-1));
energy_orig=n_fb_orig+weight_len.*tot_len_orig;

n_fb=n_fb_orig;
tot_len=tot_len_orig;
energy=energy_orig;

out=zeros(niter,1);

for i=1:niter
        
        out(i)=energy;
        i1=i_ch(floor(rand.*n_ch)+1);
        level_prev=levels(i1);
        levels(i1)=floor((nlevels).*rand)+1;
        
        n_fb_try=sum(levels(a)<=levels(b));
        tot_len_try=sum((levels(a)>levels(b)).*(levels(a)-levels(b)-1));
        energy_try=n_fb_try+weight_len.*tot_len_try;
        
        if rand>exp(-beta.*(energy_try-energy));
                levels(i1)=level_prev;
        else
                n_fb=n_fb_try;
                tot_len=tot_len_try;
                energy=energy_try;
        end;
end;

return;
