function OUTPUT_TABLE = CalcZ(root_github_dir, Aggregation_method,expt_type)
%CALCZ Calculate Z'-factor from HitTable
%   CalcZ(['per_image'|'per_well'], ['CP_CPA'|'ilastik_CP'])

close all
disp(' ') %% Just to separate the disp output
%%%%%%%%%%
%% SET THESE FILES TO THE OUTPUT FILES EXPORTED FROM DATABASE PER_IMAGE TABLES

%% We are not using 'per_image' aggregation - commenting out
% file_CPCPA_per_image = fullfile(root_github_dir,...
%     'CP_CPA/SeedingDensities/CP_CPA_IllumCorrInput_GHAD2-D6-20x/Hit_Table_Per_Image_60rules_2015_09_03.csv');
% 
% % file_ilastik_CP_per_image = fullfile(root_github_dir,'ilastik_CP/SeedingDensities/ilastik_CP_IllumCorrInput_GHAD2-D6-20x/ilastik_CP_counts_Per_Image.csv');
% file_ilastik_CP_per_image = fullfile(root_github_dir,...
%     'ilastik_CP/SeedingDensities/ilastik_CP_IllumCorrInput_GHAD2-D6-20x/Hit_Table_ilastik_CP_Per_Image_60Rules.csv');

file_CPCPA_per_well = fullfile(root_github_dir,...
    'CP_CPA/SeedingDensities/CP_CPA_IllumCorrInput_GHAD2-D6-20x/Hit_Table_Per_PlateWell_60rules_2015_09_03.csv');

% file_ilastik_CP_per_well = fullfile(root_github_dir,'ilastik_CP/SeedingDensities/ilastik_CP_IllumCorrInput_GHAD2-D6-20x/ilastik_CP_counts_Per_Well.csv');
file_ilastik_CP_per_well = fullfile(root_github_dir,...
    'ilastik_CP/SeedingDensities/ilastik_CP_IllumCorrInput_GHAD2-D6-20x/Hit_Table_ilastik_CP_Per_Well_60Rules.csv');

MAX_Y_AXIS = 1100;  %% USER SET - to normalize boxplot scales. Found empirically for per_well figures.

%%%%%%%%%%

%% CHECK INPUTS
narginchk(0,3)
if nargin < 3
    disp('Not enough inputs!  Using defaults...')
    Agg_method = 'per_image';
    expt_type = 'CP_CPA';
elseif nargin == 3
    validStr = validatestring(Aggregation_method,{'per_image','per_well'});
    Agg_method = validStr; %% Fix input if only partial match
    validStr = validatestring(expt_type,{'CP_CPA','ilastik_CP'});
    expt_type = validStr; %% Fix input if only partial match
end
disp(['Aggregation method = ' Agg_method])
disp(['Expt Type = ' expt_type])

%% GET FILENAME
switch Agg_method
    case 'per_image'
        switch expt_type
            case 'CP_CPA'
                filename = file_CPCPA_per_image;
            case 'ilastik_CP'
                filename = file_ilastik_CP_per_image;
        end
    case 'per_well'
        switch expt_type
            case 'CP_CPA'
                filename = file_CPCPA_per_well;
            case 'ilastik_CP'
                filename = file_ilastik_CP_per_well;
        end

end
disp(['Importing from ' filename])

%% Need Image_Metadata_Well, HepCellCount
% OUTPUT_TABLE = import_CSV_Zprime(Agg_method, expt_type);
OUTPUT_TABLE = readtable(filename);  

Image_Metadata_Well = OUTPUT_TABLE.Image_Metadata_Well;
HepCellCount = OUTPUT_TABLE.HepCellCount;

%% Assume X01 are the 4000 Heps and the X23 are the 9500 Heps
x_loc = [4000;9500];
% x_loc_cellstr = cellstr(num2str(x_loc)); %% Numbers for x-label
x_loc_cellstr = {'Low','High'}; %% Manually label

%% Find Control Wells (assumed to be cols '01' and '23')
charWell = char(Image_Metadata_Well);  %% COnvert to char.  OK, since they are all 3 char long
wellRow = charWell(:,2:end); %% e.g. the '03' in 'A03'
wellRowCell = cellstr(wellRow); %% Convert back to cell array, because 'strcmp' expect it
NegCtrls = strcmp(wellRowCell,'01'); %% Neg Ctrls are in col '01'
PosCtrls = strcmp(wellRowCell,'23');%% Pos Ctrls are in col '23'

HepCountsNegCtrls = HepCellCount(NegCtrls);
HepCountsPosCtrls = HepCellCount(PosCtrls);

%%%%% TEMP!! REMOVE OUTLIER
% HepCountsNegCtrls = HepCountsNegCtrls(2:end);


%% ALL [Pos|Neg]Ctrls together 
[Z_prime,meanPos,meanNeg,stdPos,stdNeg] = calc_z_prime(HepCountsNegCtrls,HepCountsPosCtrls);

%% FIGS
figure(1)
%% To get boxplot to plot groups, we need to get a grouping variable
%% And so we add concatenate the data
boxdata = [HepCountsNegCtrls; HepCountsPosCtrls];
boxgroup = [zeros(length(HepCountsNegCtrls),1); ones(length(HepCountsPosCtrls),1)];
boxplot(boxdata,boxgroup,'labels',x_loc_cellstr);
ax = axis;
if strcmp(Agg_method,'per_well')
    ax(4) = MAX_Y_AXIS;
    set(gca,'YTick',[0 200 400 600 800 1000])
end
axis([ax(1:2) 0 ax(4)]) %% Make y-axis start at zero
set(gca,'FontSize',24,'FontWeight','bold')
if strcmp(Agg_method,'per_well'), agg_text = 'per well', else agg_text = 'per image'; end
xlabel(sprintf('Hepatocytes %s',agg_text),'FontSize',32)
ylabel('Hepatocyte Counts','FontSize',32)
% title(['Hep counts ' Agg_method ';  Z'' = ' num2str(Z_prime)])
box off
axis square

% %% BARPLOT (DECIDED BOXPLOT BETTER, SO COMMENT THIS OUT FOR NOW)
% figure(2)
% bar(x_loc,[meanNeg,meanPos])
% hold on
% er = errorbar(x_loc,[meanNeg,meanPos],[stdNeg,stdPos]);
% er.Color = [0,0,0];
% er.LineStyle = 'none';
% hold off
% ylabel('Mean Hep Counts')
% set(gca,'XTickLabel',{'Low','High'}) %% Manually re-label x-axis
% title(['Mean Hep counts ' Agg_method ';  Z'' = ' num2str(Z_prime)])

% [HepCountsNegCtrls,HepCountsPosCtrls]

disp(['Z prime (Hep counts) = ' num2str(Z_prime)])

end

function [Z_prime,meanPos,meanNeg,stdPos,stdNeg] = calc_z_prime(NegCtrls,PosCtrls)
% 
USE_TRIMMEAN = 0;  %% Try trimmean, if ~= 0
TRIMMEAN_PERCENT = 5;

if USE_TRIMMEAN
    meanPos = trimmean(PosCtrls,TRIMMEAN_PERCENT);
    meanNeg = trimmean(NegCtrls,TRIMMEAN_PERCENT);
else
    meanPos = mean(PosCtrls);
    meanNeg = mean(NegCtrls);
end

stdPos = std(PosCtrls);
stdNeg = std(NegCtrls);

numer = 3.0 * (stdPos + stdNeg);
denom = abs(meanPos - meanNeg);

Z_prime = 1 - numer./denom;

end

% function OUTPUT_TABLE = import_CSV_Zprime(filename)
% %% Aggregation method should be ['per_image'|'per_well']
% %% expt_type should be ['CP_CPA'|'ilastik_CP']
% 
% 
% 
% %% USE NEW MATLAB TABLE OBJECT to extract relevant cols
% % VarNames = TBL.Properties.VariableNames;
% 
% %% This will produce a warning "Warning: Variable names were modified to make them valid MATLAB identifiers."
% %% if there are any spaces or other odd characters in the header.  It's ok, as long as the next steps looks for the 
% %% altered header names.
% OUTPUT_TABLE = readtable(filename);  
% 
% %% FOR CP Export* modules' output
% %% VarNames = TBL.Properties.VariableNames;
% %% The 
% % cols_to_keep = strncmp('Count_',TBL.Properties.VariableNames,6) | ...
% %     strncmp('FileName_',TBL.Properties.VariableNames,9) | ...
% %     strncmp('Metadata_',TBL.Properties.VariableNames,9) | ...
% %     strncmp('ObjectsFileName_',TBL.Properties.VariableNames,16) | ...
% %     strncmp('ObjectsPathName_',TBL.Properties.VariableNames,16) | ...
% %     strncmp('Overlap_',TBL.Properties.VariableNames,8) | ...
% %     strncmp('PathName_',TBL.Properties.VariableNames,9);
% % OUTPUT_TABLE = TBL(:,cols_to_keep);
% 
% %% For CPA counts output
% % cols_to_keep = strncmp('HepCellCount',TBL.Properties.VariableNames,length('HepCellCount')) | ...
% %     strncmp('Image_Metadata_Well',TBL.Properties.VariableNames,length('Image_Metadata_Well'));
% % OUTPUT_TABLE = OUTPUT_TABLE(:,cols_to_keep);
% 
% end