function [Q,QQ]=metastab(X,ci)
%   Compute the metastability according to macro-clustering ci
%       [Q,QQ]=metastab(X,ci)

numMacro=max(ci);

for i=1:numMacro,
    id{i}=find(ci==i);
end

QQ = zeros(numMacro,numMacro);
for i=1:numMacro,
    for j=i:numMacro,
        QQ(i,j)=sum(sum(X(id{i},id{j})));
        QQ(j,i)=sum(sum(X(id{j},id{i})));
    end
end

%QQ = QQ + tril(QQ',-1);

D = sum(QQ,2);
Q = sum(diag(diag(1./D)*QQ));