function [fassign,T, assignMap] = fuzzy_assign(X, gradFlowEQ, nodeInfo, opt)
% Fuzzy assignment of macrostates
%	[fassign, T, assignMap] = fuzzy_assign(X, gradFlowEQ, nodeInfo, opt)
%

if nargin<4,
	opt.fname = 'fuzzy_assign.data';
end

numMicro = max(size(X));
numNode = length(nodeInfo);

MacroStates = find(diag(gradFlowEQ)==1);
numMacro = length(MacroStates);

% Compute the implied equilibrium distribution from each node. 
distEQ = gradFlowEQ(MacroStates,:);

fass = zeros(numMicro,numMacro);

for i=numNode:-1:1,
	fass(nodeInfo{i}.set,:)=repmat(distEQ(:,i)',[length(nodeInfo{i}.set),1]);
end

% compute zero rows as random assignment
id = find(sum(fass,2)==0);
fass(id,:) = ones(length(id),numMacro)./numMacro;

mass = max(fass,[],2);
id_fuzzy = find(mass < 1);


% map the original assignment to the new
assignMap = zeros(numMicro,1);
assignMap(id_fuzzy) = [numMacro+1:1:numMacro+length(id_fuzzy)];

% Construct new transition count matrix from X
T = speye(numMacro+length(id_fuzzy));
for i=1:numMacro,
	id_nonfuzzy{i}=find(fass(:,i)==1);
	assignMap(id_nonfuzzy{i})=i; % assignment mapping
end

for i=1:numMacro,
	for j=i:numMacro,
		T(i,j) = sum(sum(X(id_nonfuzzy{i},id_nonfuzzy{j})));
	end
	for j=1:length(id_fuzzy),
		T(i,j+numMacro) = sum(X(id_nonfuzzy{i},id_fuzzy(j)));
	end
end
T = T + triu(T,1)';

T(numMacro+1:numMacro+length(id_fuzzy),numMacro+1:numMacro+length(id_fuzzy)) = X(id_fuzzy,id_fuzzy);

% rearrange the fuzzy assignment
fassign = [eye(numMacro);fass(id_fuzzy,:)];
% fassign = repmat(1./sum((fassign>0),2),[1,numMacro]).*(fassign>0);
% T2 = T^2;
% T4 = T2^2;
% Tm = full(T4(numMacro+1:size(T,1),1:numMacro));
% dm = sum(Tm,2);
% fassign = [eye(numMacro);repmat(1./dm,[1,numMacro]).*Tm];

% write to the fuzzy assign files
fname = opt.fname; % 'fuzzy_assign.data';
[F_rowsort, id_rowsort]=sort(fassign,2,'descend');
fid = fopen(fname,'w');
fprintf(fid,'%d\t%d\n',size(fassign,1),numMacro);
for i=1:size(fassign,1),
	if i<=numMacro,
		fprintf(fid,'%d\t%d\t%d\t%d\n',i,1,i,1);
	else
		id = find(F_rowsort(i,:)>1e-4);
		fprintf(fid,'%d\t%d',i,length(id));
		fprintf(fid,'\t%d\t%0.4f',[id_rowsort(i,id);F_rowsort(i,id)]);
		fprintf(fid,'\n');
	end,
end

fclose(fid);

