%% Mapper with super level sets to construct the macrostates
%
%  	1) Run Mapper with super level sets to construct the graph: mapper_super.m (mapper_yy)
%	2) Analyze the gradient flow on the mapper graph: flowGrad.m (mapper_yy)
%	3) Optimal assignment of microstates by diffusions according to transition count matrix: assign_optim.m (mapper_yy)

% Yuan Yao, 4/22/09
% Stanford University

%% May need Statistics (commercial), matlab_bgl (free) toolboxes
%

if ~exist('mapper_yy','dir'),
	addpath ../mapper_yy
end

%load ../anlmatrix/alanine_dipeptide_minmax_clustering_k500_seed0.mat cluster_size orbit
load T10000_s1.mat T
X = (T+T')/2;
cluster_size = full(sum(T,2));
%filter = full(sum(X,2));
filter = cluster_size(:,1);

%randSurf = 0.1;

%%%%%% Default choice of levels 
% levels = [100, 1000, 3000 5000];
% clusters = [2, 4, 6, 6];

%%%%%% Choose Levels according to percentage of samples
plevels = [0.1:0.2:0.9];
fcdf = hcdf(filter);
for i=1:length(plevels), levels(i) = sum(fcdf<=plevels(i)); end
%clusters = [2 2 3 4 4 4 4 4 6];
clusters = zeros(length(plevels),1);

%%%%%% Choose levels according to number of microstates
%levels = [100:100:2500];

numLevels = length(levels);
[b,g]=apBetti(X,filter,levels);
for i=1:numLevels,
	clusters(i) = b{i}(1);
end
superLevelSet = superLvlSet(filter,levels,clusters);

opt.clustSize = 0.01; %% Keep those clusters larger than 1% of the level size 

% 1) Call mapper_super to construct the connectivity graph
[adja, nodeInfo, levelIdx] = mapper_super(X, superLevelSet); % opt

% 2) Analyze the gradient flow on the Mapper graph
cpt = flowGrad(adja,levelIdx,nodeInfo); 

label{1} = sprintf(' Number of Levels : %d', length(superLevelSet));
disp('Mapper done.');

fname = sprintf('RNA_k%d_L%d_M%d',length(X),numLevels,length(cpt.localMax));

save(sprintf('%s.mat',fname));

% write gradient flow to .dot file
DotFileName = sprintf('%s_flow',fname);
writeFlowGraph(sprintf('%s.dot',DotFileName), cpt.gradFlow, nodeInfo ,superLevelSet,levelIdx); %% a function in mapper toolbox
disp(sprintf('%s.dot File done.',DotFileName));

% convert .dot to .jpg file
system(sprintf('dot -Tjpg %s.dot > %s.jpg', DotFileName, DotFileName)); %% need Graphviz
disp(sprintf('%s.jpg File done.',DotFileName));
	
% write node information to a .txt file
%	1st row: number of nodes
%	each row followed is a list of microstates associated with a node
NodeFileName = sprintf('%s_node.txt',fname);
fid = fopen(NodeFileName,'w');
fprintf(fid,'%d \n',length(nodeInfo)); % 1st row: number of node
for i=1:length(nodeInfo),
	fprintf(fid,'%d ',nodeInfo{i}.set); % each row is a list of microstate ids
	fprintf(fid,' \n');
end
fclose(fid);

% assignFile = sprintf('%s_fuzzy_assign.txt',fname);
% opt.fname = assignFile;
% [fassign,T,assignMap]=fuzzy_assign(X,cpt.equilibrium, nodeInfo,opt);

% 3) Get optimal assignment of microstates by Diffusion according to X, the transition count matrix 
[ci,csize,fassign,T,Qmax,id_fuzzy]=assign_optim(X,cpt.equilibrium,nodeInfo,opt);

% T = full(T);
% transitFile = sprintf('%s_fuzzy_transit.txt',fname);
% save(transitFile,'T','-ascii');
% outFile = sprintf('ala_sa3_super_optimal_assign_k%d_L%d_%d.txt',length(X),numLevels,length(T));
% logFile = sprintf('ala_sa3_super_k%d_L%d_%d.out',length(X),numLevels,length(T));
% system(sprintf('../sa/src/run_SA3 -t %s -a %s -o %s > %s &',transitFile,assignFile,outFile,logFile));
% imageFile = sprintf('ala_super_M%d_image_k%d',size(fassign,2),length(X));

disp('Done.');
% disp('  Use "type(logFile)" to check SA Run in background.');
% disp('  After that, use "imshowMacro" to check the pictures of macrostates');

save(sprintf('%s.mat',fname));
