clear all
close all

addpath(genpath('/data'))

try
    load('biomodels_modeldims.mat')
    imodel = length(nx) + 1;
catch
    filelist = dir('curated');
    
    model = TranslateSBML('./curated/BIOMD0000000255.xml',0,0);
    modelUniProtIDs = unique(cellfun(@(x) extractUniProtIdentifier(x),{model.species.annotation}));
    modelUniProtIDs(modelUniProtIDs==0) = [];
    
    imodel = 1;
    for ifile = 1:length(filelist)
        if(filelist(ifile).bytes>0)
            try
                model = TranslateSBML(fullfile(filelist(ifile).folder,filelist(ifile).name),0,0);
                
                nr(imodel) = length(model.reaction);
                nx(imodel) = sum(not([model.species.boundaryCondition]));
                UniProtIDs{imodel} = unique(cellfun(@(x) extractUniProtIdentifier(x),{model.species.annotation}));
                UniProtIDs{imodel}(UniProtIDs{imodel}==0) = [];
                name{imodel} = model.name;
                overlap(imodel) = sum(ismember(UniProtIDs{imodel},modelUniProtIDs));
                missing(imodel) = sum(not(ismember(UniProtIDs{imodel},modelUniProtIDs)));
                imodel = imodel + 1
            catch err
                disp(err.message)
            end
        end
    end
    % examples
    
    % Zangh 2014
    nr(imodel) = 60;
    nx(imodel) = 59;
    cancer(imodel) = 2;
    
    imodel = imodel+1;
    
    % Chen 2009
    
    model = TranslateSBML('./curated/BIOMD0000000255.xml',0,0);
    nr(imodel) = length(model.reaction);
    nx(imodel) = sum(not([model.species.boundaryCondition]));
    
    cancer(imodel) = 2;
    
    imodel = imodel+1;
    
    % Hass 2017
    nr(imodel) = 113;
    nx(imodel) = 61;
    
    cancer(imodel) = 2;
    
    imodel = imodel+1;
    
    % Kholodenko 1999
    
    nr(imodel) = 20;
    nx(imodel) = 25;
    
    cancer(imodel) = 2;
    
    imodel = imodel+1;
    
    model = TranslateSBML('./../supplement_code/NBT_Signalling_ERBB_RAS_AKT.xml',0,0);
    nr(imodel) = length(model.reaction);
    nx(imodel) = sum(not([model.species.boundaryCondition]));
    
    cancer(imodel) = 3;
    
    imodel = imodel+1;
    
    save('biomodels_modeldims','nr','nx','cancer','name','missing','overlap','UniProtIDs');
end

% Bouhaddou 2018
nr(imodel) = 158+34+190+11+24+66+104+87;
nx(imodel) = 774;

cancer(imodel) = 2;

imodel = imodel+1;



figure
g = gramm('x',nr,'y',nx,'color',cancer)
g.geom_point();
g.axe_property('YScale','log','XScale','log','XLim',[10^(-0.5),10^4.5],'YLim',[10^(-0.5),10^4.5],'FontSize',12,'LineWidth',1,'DataAspectRatio',[1,1,1],'XTick',logspace(0,4,3),'YTick',logspace(0,4,3))
g.no_legend();
g.set_names('y','number of species','x','number of reactions')
g.draw();

makeFigure(0,0,2,150,130,'figure_Biomodels',true)
%%
figure
x = {'Bouhaddou','Bouhaddou',...
    'proposed model','proposed model',...
    'Zhang2014','Zhang2014',...
    'Chen2009','Chen2009',...
    'Hess2017','Hess2017',...
    'Kholodenko1999','Kholodenko1999',...
    'Capuani2015','Capuani2015'}
    
y = [41,100,...
     108,0,...
     14,5,...
     17,1,...
     9,6,...
     6,0,...
     4,0]
     
color = [kron(ones(1,7),[1,0])];
g = gramm('x',x,'y',y,'color',color);
g.geom_bar('stacked',true);
g.set_order_options('x',{'Bouhaddou','proposed model','Zhang2014','Chen2009','Hess2017','Kholodenko1999','Capuani2015'},'color',[1,0]);
g.axe_property('XTickLabelRotation',90);
g.set_names('y','number of genes/proteins','x','model')
g.no_legend;
g.draw;

makeFigure(0,0,2,70,135,'figure_models_genes',true)
     


%
% model = TranslateSBML('Model_ERBB.xml')
% compartments = {model.species.compartment};
% compartments_noBC = compartments(not([model.species.boundaryCondition]));
% sum(strcmp('compartment_2',compartments_noBC))
% sum(strcmp('compartment_3',compartments_noBC))
% sum(strcmp('compartment_5',compartments_noBC))
% sum(strcmp('compartment_5',compartments_noBC))


function identifier = extractUniProtIdentifier(xml)
    tok = regexp(xml,'identifiers\.org\/uniprot\/P([0-9]*)','tokens');
    if(~isempty(tok))
        identifier = str2double(tok{1}{1});
    else
        identifier = 0;
    end
end