clear all
prep_figure
close all

%%
plotConcFun(sim_fit,data_fit,conc,drugs,'corr','fit - corr');
makeFigure(0,0,2,180,120,'mechanistic_corr_fit',true)
plotConcFun(sim_val,data_val,conc,drugs,'corr','val - corr');
makeFigure(0,0,2,180,120,'mechanistic_corr_val',true)
plotConcFun(sim_test,data_test,conc,drugs,'corr','test - corr');
makeFigure(0,0,2,180,120,'mechanistic_corr_test',true)
%%
[xx_mechanistic_fit,yy_mechanistic_fit] = plotConcFun(sim_fit,data_fit,conc,drugs,'mse','fit - rmse');
makeFigure(0,0,2,180,120,'mechanistic_rmse_fit',true)
[xx_mechanistic_val,yy_mechanistic_val] = plotConcFun(sim_val,data_val,conc,drugs,'mse','val - rmse');
makeFigure(0,0,2,180,120,'mechanistic_rmse_val',true)
[xx_mechanistic_test,yy_mechanistic_test] = plotConcFun(sim_test,data_test,conc,drugs,'mse','test - rmse');
makeFigure(0,0,2,180,120,'mechanistic_rmse_test',true)

%%
tmp = permute(median(data_fit,5),[4,1,2,3]);
nanvar(tmp(:,:),[],2)

tmp = permute(median(data_val,5),[4,1,2,3]);
nanvar(tmp(:,:),[],2)

tmp = permute(median(data_test,5),[4,1,2,3]);
nanvar(tmp(:,:),[],2)

function [xx,yy] = plotConcFun(sim,data,conc,drugs,mode,title)
    
    switch(mode)
        case 'corr'
            fun = @(x,y) getcorrcoeff(x,y);
            ylim = [-1,1];
        case 'mse'
            fun = @(x,y) getmse(x,y);
            ylim = [0,0.5];
    end
    
    drugs{end+1} = 'all';
    for ic = 1:size(data,3)
        for icv = 1:5
            for id = 1:size(data,2)
                x = nanmedian(sim(:,id,ic,icv,:),5);
                y = nanmedian(data(:,id,ic,icv,:),5);
                z = fun(x(~isnan(y)),y(~isnan(y)));
                
                xx(id,ic,icv) = ic;
                yy(id,ic,icv) = z;
                cc(id,ic,icv) = id;
            end
            
            x = nanmedian(sim(:,:,ic,icv,:),5);
            y = nanmedian(data(:,:,ic,icv,:),5);
            z = fun(x(~isnan(y)),y(~isnan(y)));
            
            id=id+1;
            
            xx(id,ic,icv) = ic;
            yy(id,ic,icv) = z;
            cc(id,ic,icv) = id;
        end
    end
    figure
    g = gramm('x',conc(xx(:))*1000,'y',yy(:),'color',drugs(cc(:)));
    if(numel(unique(xx(:)))>1)
        g.stat_summary('type','ci','geom','lines');
        g.axe_property('XScale','log','YLim',ylim,'XLim',[2,1e4]);
        switch(mode)
            case 'corr'
                g.set_names('x','drug concentration [nM]','y','Pearson correlation coefficient','color','drug');
            case 'mse' 
                g.set_names('x','drug concentration [nM]','y','root-mean-square error','color','drug');
        end
    else
        g.stat_boxplot()
        g.set_names('x','drug concentration','y','correlation coefficient','color','drug');
    end

    g.set_title(title);
    g.draw;
end

function z = getcorrcoeff(x,y)
    [C,P] = corrcoef(x(~isnan(y)),y(~isnan(y)));
    z=C(1,2);
end

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