clear all
prep_figure
close all

%%
index = 1;
crossVar = [];
selTissues = {};
selDrugs = {};
modes = {};
values = [];


for mode = {'rmse','corr'}
    for icv = 1:5
        for itissue = 1:length(tissues)
            x = nanmedian(data_fit(tissue==itissue,:,:,icv,:),5);
            y = nanmedian(sim_fit(tissue==itissue,:,:,icv,:),5);
            switch(mode{1})
                case 'corr'
                    values(index) = getcorrcoeff(x,y);
                case 'rmse'
                    values(index) = getrmse(x,y);
            end
            crossVar(index) = icv;
            selTissues{index} = tissues{itissue};
            selDrugs{index} = 'all of the conditions';
            modes{index} = mode{1};
            index = index + 1;
        end
        for idrug = 1:7
            x = nanmedian(data_fit(:,idrug,:,icv,:),5);
            y = nanmedian(sim_fit(:,idrug,:,icv,:),5);
            switch(mode{1})
                case 'corr'
                    values(index) = getcorrcoeff(x,y);
                case 'rmse'
                    values(index) = getrmse(x,y);
            end
            crossVar(index) = icv;
            selTissues{index} = 'all of the conditions';
            selDrugs{index} = drugs{idrug};
            modes{index} = mode{1};
            index = index + 1;
        end
        x = nanmedian(data_fit(:,:,:,icv,:),5);
        y = nanmedian(sim_fit(:,:,:,icv,:),5);
        switch(mode{1})
            case 'corr'
                values(index) = getcorrcoeff(x,y);
            case 'rmse'
                values(index) = getrmse(x,y);
        end
        crossVar(index) = icv;
        selTissues{index} = 'all of the conditions';
        selDrugs{index} = 'all of the conditions';
        modes{index} = mode{1};
        index = index + 1;
    end
end

%%
figure
for iplot = 1:2
    switch(iplot)
        case 1
            xx = selTissues;
            subset = strcmp(selDrugs,'all of the conditions');
        case 2
            xx = selDrugs;
            subset = strcmp(selTissues,'all of the conditions');
    end
    g(iplot) = gramm('x',xx,'y',values,'subset',subset)
    g(iplot).facet_grid(modes,[],'scale','free_y')
    g(iplot).stat_summary('geom','bar')
    g(iplot).stat_summary('geom','errorbar')
    g(iplot).axe_property('XTickLabelRotation',90)
end

g.draw()
makeFigure(0,0,2,170,140,'mechanistic_corr_rmse_tissue_drugs',true)

%%

for mode = {'rmse','corr'}
    mode{1}
    index = 1;
    pvals = [];
    ref = values(all([strcmp(selDrugs,'all of the conditions'),
                           strcmp(selTissues,'all of the conditions'),
                           strcmp(modes,mode{1})]));
    mean(ref)
    std(ref)
    for idrug = 1:length(drugs)
        [~,p] = ttest(ref,values(all([strcmp(selDrugs,drugs{idrug}),
                           strcmp(selTissues,'all of the conditions'),
                           strcmp(modes,mode{1})])));
        pvals(index) = p;
        index = index + 1;
    end
    for itissue = 1:length(tissues)
        [~,p] = ttest(ref,values(all([strcmp(selDrugs,'all of the conditions'),
                           strcmp(selTissues,tissues{itissue}),
                           strcmp(modes,mode{1})])));
        pvals(index) = p;
        index = index + 1;
    end
    % BH correction
    pvals = pvals*length(pvals);

    drugs(pvals(1:7)<0.05)
    tissues(pvals(8:end)<0.05)
    pvals([3,4,11])
end


        
function z = getcorrcoeff(x,y)
    z = corr(x(~isnan(y)),y(~isnan(y)),'type','Pearson');
end

function z = getrmse(x,y)
    z = sqrt(nanmean((x(:)-y(:)).^2));
end