load('Rall.mat')
load('rnames.mat')
load('pnames.mat')


FDR = 0.05;

par = [];
for icv = 1:5
    [~,msidx] = sort(arrayfun(@(x) x.fval(end),R(:,icv)),'ascend');
    pars = arrayfun(@(x)x.par(:),R(:,icv),'UniformOutput',false);
    par = [par,horzcat(pars{msidx(1:5)})];
end

%%
Kmat = [];
K_name = {};
for idx = find(cellfun(@(x) ~isempty(strfind(x,'|k0}')),p_name));
    jdx = find(strcmp(strrep(p_name(idx),'k0','kd'),p_name));
    if(~isempty(jdx))
        Kmat = [Kmat;[idx,jdx]];
        K_name = [K_name,p_name{idx}];
    else
%         Kmat = [Kmat;[idx,0]];
%         K_name = [K_name,p_name{idx}];
    end
    
end

Kval = par(Kmat(:,1),:);
Kval(Kmat(:,2)>0,:) = Kval(Kmat(:,2)>0,:)-par(Kmat(Kmat(:,2)>0,2),:);

[h,p] = ttest(Kval',0);
[cor_p, h]=bonf_holm(p,FDR);

k0s = find(strcmp(strrep(p_name(idx),'k0','kd'),p_name));


%%
xnames=cellfun(@(x) getReactionName(x),K_name','UniformOutput',false);
xnames=regexprep(xnames,'-00[0-9]+','');
xnames=regexprep(xnames,' \[[\w ]*\]','');
xnames=regexprep(xnames,'P\[[\w;]*\]-','p');
figure
npar = length(xnames);
cc = kron(kron(1:5,ones(1,5)),ones(npar,1));
xx = kron(kron(ones(1,5),ones(1,5)),(1:npar)');
subset = kron(kron(ones(1,5),ones(1,5)),cor_p<=FDR);
yy = Kval(:,:);



clear g
g(1,1) = gramm('x',yy)
g(1,1).stat_density()
g(1,1).axe_property('XLim',[-3,3],'XTick',[-3:3],'TickLabelInterpreter','tex');

g(2,1) = gramm('x',xnames(xx(:)),'y',yy(:),'subset',logical(subset));
g(2,1).stat_boxplot;
g(2,1).coord_flip;
g(2,1).axe_property('YLim',[-3,3],'YTick',[-3:3],'TickLabelInterpreter','tex');
g(2,1).set_names('x','Reaction description','y','Equilibrium Constant')
% g(2,1).set_order_options('x',xnames)

g.draw();
% g.facet_axes_handles(1).YTickLabel = arrayfun(@(x) ['10^{' num2str(x) '}'],g.facet_axes_handles(1).YTick,'UniformOutput',false)

makeFigure(0,0,2,400,250,'figure_parameter_estimates',true)

%%


    
function name = getReactionName(parname)
    persistent names ids
    
    if(isempty(names) || isempty(ids))
        load('rnames.mat')
        names = R_name(:,2);
        ids = R_name(:,1);
    end
    
    name = 'no reaction found';
    tok = regexp(parname,'\{([0-9]*)\|','tokens');
    if(~isempty(tok))
        nameidx = find(strcmp(tok{1},ids));
        name = names{nameidx};
    end
end
    
    
