% Syntax: [levels, n_fb, tot_len]=hierarchy_levels(s1,nlevels,itop,ibot,levels,beta, niter, 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+2)
% 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
function [levels, n_fb, tot_len]=hierarchy_levels(s1,nlevels,itop,ibot,beta, niter,levels);

weight_len=0.1;

s1=sign(abs(s1-diag(diag(s1))));
n_edges=full(sum(sum(s1)));
len1=length(s1);
ap=sign(sum(s1+s1'));
a_changeable=ap;

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

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

if (nargin < 7) 
  levels=sparse(len1,1);
  ip=find(a_changeable);
  levels(ip)=floor((nlevels).*rand(size(ip)))+2;
end;


if isempty(itop)==0
    levels(itop)=nlevels+2;
    ap(itop)=1;
    a_changeable(itop)=0;
end;

if isempty(ibot)==0
    levels(ibot)=1;
    ap(ibot)=1;
    a_changeable(ibot)=0;
end;

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

[a,b,c]=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;

for i=1:niter;
    i1=i_ch(floor(rand.*n_ch)+1);
    level_prev=levels(i1);
    levels(i1)=floor((nlevels).*rand)+2;
    
    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;