function [GOnum_vec,MF_cell,graphs,idx_in_graphs,GOnum_to_index,GOidx_to_MFidx,...
    index_bad,index_good] = readsubgraphs_oct05(filename,num_graphs) 
% Reads Borya's file and records it in a comfortable matlab seeting
% Input:
% filename (character):  name of data file
% num_graphs (integer):  number of graphs in file 'filename' to read (ascending order)
% Output:
%
% GOnum_vec (real vector maxidx*1).
%                Vector of numbers in GO notation (e.g. 3824 for GO:0003824)
%                according to their order of apperance in file 'filename'
% MF_cell (cell of length num_graphs)
%                It lists the name of molecular function (MF) for each
%                graph (e.g. the first one is 12asA_1).
% graphs        (cell of length num_graphs)
%                For any given MF (molecular function) according to order in
%                file, it provides sparse matrix of graph representing that function
% idx_in_graphs (cell of length num_graphs)
%                For any given MF (molecular function) according to order in
%                file, it lists the indices of GO notation according to 
%                order in file.
% GOnum_to_index (sparse array 100000*1).
%                For each number of GO notation (e.g. 3824 for GO:0003824)
%                indexing the array, it provides the index according to the
%                order in file 'filename'% maxidx (integer). 
%                The highest number of GO notation involved, which is the largest 
%                size of graphs
% GOidx_to_MFidx (a cell of length maxidx).
%                For each index of GO notation (according to order in file) 
%                it provides all indices of MF (molecular functions) 
%                associated with it.
% index_bad (coloumn vector of length = number of empty MF graphs)
%                vector of indices (according to location in files) of all 
%                MF functions (graphs) with no GO notation (empty graphs).
% index_good (coloumn vector of length = num_graphs - # of empty MF graphs)
%                vector of indices (according to location in files) of all 
%                MF functions (graphs) with nonempty graphs.
% Remarks:
% the inner variable graphs is a cell array of size num_graphs 
% (or num_graphs + 1, if num_graphs < max number)
% of all graphs including empty ones (if num_graphs < max number, then 
% the last graph is empty.
% idx_in_graphs (cell array of indices used in each graph of length num_graphs)

if nargin < 2
    %num_graphs = 3306;
    %num_graphs = 3032;
    num_graphs = 3139;
end
fd = fopen(filename);
MF_idx = 0; %index of molecular function (MF, graph) following order in file 'filename'
GOnum_to_index = sparse(100000,1); %sparse vector to map numbers apperaing 
%in GO notation to indices according to the order of notation in file
                     
MF_cell = cell(num_graphs,1);

GO_idx = 0; %index to GO notation according to order in file 'filename'
                
if(fd == -1)
    error('Error opening file!');
    return;
end

s = fgetl(fd);
while(ischar(s)) & MF_idx<=num_graphs
    if(length(s) > 0)
        switch(s(1))
            case 'D' % a line announcing a new graph
                MF_idx = MF_idx+1;
                graphs{MF_idx} = sparse(1,1);
                idx_in_graphs{MF_idx} = 0; %dummy choise in order not to stop in 
                                         %if length(setdiff(idx_in_graphs{MF_idx},GOnum1))...
                internal_idx = 0;
                if(mod(MF_idx,1000) == 0)
                    disp('read graphs');
                    disp(MF_idx);
                end
                prot = '\w*_\w*';
                m = regexp(s, prot, 'match');
                MF_cell{MF_idx} = m{1};
            case 'G' % a line of two GO notation and weight (e.g. GO:0004888 GO:0004872	0.04761905)
                %process line.
                [start ends extents match tokens name] = regexp(s,'((\.|[1-9])[0-9]*)');
                GOnum1 = sscanf(match{1},'%d');
                GOnum2 = sscanf(match{2},'%d');
                %weight = sscanf(match{3},'%f');
                weight = 1;
                %We need to create a new index to go from index in file 'filename' to
                %the index we will use.
                %  break this off into another function?

                idx1 = GOnum_to_index(GOnum1);
                if(idx1 == 0)
                    GO_idx = GO_idx + 1;
                    idx1 = GO_idx;
                    GOnum_to_index(GOnum1) = idx1;
                    GOnum_vec(idx1) = GOnum1;
                    GOidx_to_MFidx{idx1} = MF_idx;
                elseif ~ismember(MF_idx,GOidx_to_MFidx{idx1})
                    GOidx_to_MFidx{idx1} = [GOidx_to_MFidx{idx1},MF_idx];
                end

                idx2 = GOnum_to_index(GOnum2);
                if(idx2 == 0)
                    GO_idx = GO_idx + 1;
                    idx2 = GO_idx;
                    GOnum_to_index(GOnum2) = idx2;
                    GOnum_vec(idx2) = GOnum2;
                    GOidx_to_MFidx{idx2} = MF_idx;
                elseif ~ismember(MF_idx,GOidx_to_MFidx{idx2}) 
                    GOidx_to_MFidx{idx2} = [GOidx_to_MFidx{idx2},MF_idx];
                end
                
                
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      
      
                %graphs{MF_idx}(idx1,idx2) = 1;
                %graphs{MF_idx}(idx2,idx1) = 1;
                graphs{MF_idx}(idx1,idx2) = weight;
                graphs{MF_idx}(idx2,idx1) = weight;
                if length(setdiff(idx_in_graphs{MF_idx},idx1)) == length(idx_in_graphs{MF_idx})
                    internal_idx = internal_idx + 1;
                    idx_in_graphs{MF_idx}(internal_idx) = idx1;
                end
                if length(setdiff(idx_in_graphs{MF_idx},idx2)) == length(idx_in_graphs{MF_idx})
                    internal_idx = internal_idx + 1;
                    idx_in_graphs{MF_idx}(internal_idx) = idx2;
                end
        end
    end
    s = fgetl(fd); % get next line
end

%%%%%%%% Fill out sparse matricies to same sizes and sort out empty graphs;
%if num_graphs < max number, we get an adittinal graph of index num_graphs+1
%whos matrix is zeros(maxidx,maxidx) so it will appear with the empty ones

FILL = 1;
good_MF_idx = 0;
bad_MF_idx = 0;


%bad_cell = [];
index_bad = [];
if(FILL)    
  maxidx = GO_idx;
  for i=1:MF_idx
    if sum(idx_in_graphs{i})==0
       bad_MF_idx = bad_MF_idx + 1; 
       index_bad(bad_MF_idx) = i;
    else
       good_MF_idx = good_MF_idx + 1;
       index_good(good_MF_idx) = i;
    end
       graphs{i}(maxidx,maxidx) = 0;
  end
end
if size(GOnum_vec,2)>1
    GOnum_vec = GOnum_vec';
end
%bad_cell = MF_cell(index_bad);
%good_cell = MF_cell(index_good);
%good_graphs = graphs(index_good);
%goodidx_in_graphs = idx_in_graphs(index_good);