clear all; close all;

feature = 'charge';

%cultureID = {'anchovy5','anchovy4','anchovy6'};
%culturelabels = {'control','TTX','TTX + CTZ'};
% cultureID = {'dingbat4','dingbat3'};
% culturelabels = {'control','CTZ'};

% cultureID = {'stimpy2','stimpy3','stimpy4'};
% culturelabels = {'control','TTX','TTX + CTZ'};
% dirID = 'C:\Users\mingfaifong\Documents\2013-01-05_stimpy2-stimpy3-stimpy4_patch-analysis\';
% savename = [dirID,'2013-01-05_stimpy2-stimpy3-stimpy4_', feature];

% cultureID = {'stimpy1','stimpy6','stimpy5'};
% culturelabels = {'control','TTX','TTX + CTZ'};
% dirID = 'C:\Users\mingfaifong\Documents\2013-01-04_stimpy1-stimpy6-stimpy5_patch-analysis\';
% savename = [dirID,'2013-01-05_stimpy1-stimpy6-stimpy5_', feature];

% cultureID = {'usain2','usain22493','usain22471'};
% culturelabels = {'control','CNQX','CNQX+stim'};
% dirID = 'C:\Users\mingfaifong\Documents\2013-01-30_usain2-usain22493-usain22471_patch-analysis\';
% savename = [dirID,'2013-01-30_usain2-usain22493-usain22471_', feature];

% cultureID = {'trouble1','trouble22467','trouble22468'};
% culturelabels = {'control','CNQX','CNQX+stim'};
% dirID = 'C:\Users\mingfaifong\Documents\2013-01-23_trouble1-trouble22467-trouble22468_patch-analysis\';
% savename = [dirID,'2013-01-23_trouble1-trouble22467-trouble22468_', feature];

% cultureID = {'wingnut1','wingnut22468','wingnut21332'};
% culturelabels = {'control','CNQX','CNQX+stim'};
% dirID = 'C:\Users\mingfaifong\Documents\2013-02-05_wingnut1-wingnut22468-wingnut21332_patch-analysis\';
% savename = [dirID,'2013-02-05_wingnut1-wingnut22468-wingnut21332_', feature];
% 
% cultureID = {'xinhai1','xinhai22493','xinhai22471'};
% culturelabels = {'control','CNQX','CNQX+stim'};
% dirID = 'C:\Users\mingfaifong\Documents\2013-02-10_xinhai1-xinhai22493-xinhai22471_patch-analysis\';
% savename = [dirID,'2013-02-10_xinhai1-xinhai22493-xinhai22471_', feature];

% cultureID = {'zarathustra2','zarathustra22468','zarathustra21332'};
% culturelabels = {'control','CNQX','CNQX+stim2x'};
% dirID = 'C:\Users\mingfaifong\Documents\2013-02-21_zarathustra2-zarathustra22468-zarathustra21332_patch-analysis\';
% savename = [dirID,'2013-02-21_zarathustra2-zarathustra22468-zarathustra21332_', feature];

% cultureID = {'barbie1','barbie2','barbie22468'};
% culturelabels = {'control','CNQX','CNQX+stim2x'};
% dirID = 'C:\Users\mingfaifong\Documents\2013-03-10_barbie1-barbie2-barbie22468_patch-analysis\';
% savename = [dirID,'2013-03-10_barbie1-barbie2-barbie22468_', feature];

% cultureID = {'atlas2','atlas22493'};
% culturelabels = {'control','CNQX'};
% dirID = 'C:\Users\mingfaifong\Documents\2013-03-03_atlas2-atlas22493_patch-analysis\';
% savename = [dirID,'2013-03-03_atlas2-atlas22493_', feature];

cultureID = {'dojo1','dojo22493','dojo22471'};
culturelabels = {'control','CNQX','CNQX+stim2x'};
dirID = 'C:\Users\mingfaifong\Documents\2013-03-21_dojo1-dojo22493-dojo22471_patch-analysis\';
savename = [dirID,'2013-03-21_dojo1-dojo22493-dojo22471_', feature];


% cultureID = {'juju1','jujuI5372','jujuI5388'};
% culturelabels = {'control','CNQX','CNQX+stim'};
% dirID = 'C:\Users\mingfaifong\Documents\2012-10-30_juju1-jujuI5372-jujuI5388_patch-analysis\';
% savename = [dirID,'tmp_stats\2012-10-30_juju1-jujuI5372-jujuI5388_', feature];



% cultureID = {'xeko1','xeko20629','xeko20628'};
% culturelabels = {'control','CNQX','CNQX + stimulation'};
% savename = ['C:\Documents and Settings\mingfaifong\My Documents\MiniAnalysis Data\', ...
%     '2012-05-24_xeko1-xeko20629-xeko20628_', feature];

n_cutoff = 30;
x_offset = 0;

plotcolors = {[0 0 0],[0 0.2 0.7],[0.9 0 0],[0.1 0.7 0.1],[1 0.5 0],[0.6 0 0.7]}; %black, blue, red, green, orange, purple
%plotcolors = {[0 0 0],[0.1 0.7 0.1],[1 0.5 0],[0.6 0 0.7]}; %black, green, orange, purple

% -- setting values for feature of interest -- %
switch feature
    case 'amplitude'
        f_label = 'Amplitude';
        f_unit = 'pA';
        f_datcol = 3;
        f_bounds = [5 400];
        f_cumhist = 1;
    case 'frequency'
        f_label = 'Frequency';
        f_unit = 'Hz';
        f_datcol = 19;
        f_bounds = [1/200 1000];
        f_cumhist = 0;
    case 'decay'
        f_label = 'Decay Time Constant';
        f_unit = 'ms';
        f_datcol = 5;
        f_bounds = [0.2 6]; %for AMPAergic minis
        f_cumhist = 0;
    case 'charge'
        f_label = 'Charge Transfer';
        f_unit = 'fC';
        f_datcol = 6;
        f_bounds = [1 inf];
        f_cumhist = 0;
end

f1 = figure(1);
set(f1,'Position',[1 1 1000 700])
f2 = figure(2);
set(f2,'Position',[100 200 500 400])
f3 = figure(3);
set(f3,'Position',[600 200 500 400])


%% plot by cell, display averages %%

byCell_ticklabels = {''};
disp([10,'============================================',10,f_label ' (in ' f_unit '), mean +/- SEM',10,'============================================'])

tmp_plotdata = [];
for i_cultureID = 1:length(cultureID)
    load([dirID cultureID{i_cultureID} '_alldata.mat']);
    plotdata_all = [];
    plotdata_good = [];
    plotdata_bad = [];
    clearvars byCell* -except byCell_ticklabels
    
    % -- generating data for plotting -- %
    for i_dat = 1:length(dat)
        x = dat{i_dat}(:,f_datcol);
        x = x(x>=f_bounds(1) & x<=f_bounds(2));
        tmp_plotdata = [i_dat*ones(size(x)),x]; %first column is cell number, second column is feature data
        plotdata_all = [plotdata_all; tmp_plotdata];
        if length(tmp_plotdata) >= n_cutoff
            plotdata_good = [plotdata_good; tmp_plotdata];
            byCell_good(i_dat) = 1;
        elseif isempty(tmp_plotdata)
            byCell_cellID(i_dat) = i_dat;
            byCell_ticklabels{end+1} = i_dat;
            continue
        else
            plotdata_bad = [plotdata_bad; tmp_plotdata];
            byCell_good(i_dat) = 0;
        end
        byCell_avg(i_dat) = nanmean(tmp_plotdata(:,2));
        byCell_n(i_dat) = length(tmp_plotdata(:,2));
        byCell_nanstderr(i_dat) = nanstd(tmp_plotdata(:,2))/sqrt(byCell_n(i_dat));
        byCell_cellID(i_dat) = i_dat;
        byCell_ticklabels{end+1} = i_dat;
    end

    % -- displaying cell and overall averages -- %
    disp([10,cultureID{i_cultureID},': ',culturelabels{i_cultureID},10,'--------------------------------------------'])
    for i_avg = 1:length(byCell_avg)
        disp(['cell ' num2str(i_avg) ': ' num2str(byCell_avg(i_avg)) ' +/- ' num2str(byCell_nanstderr(i_avg)) ' (n=' num2str(byCell_n(i_avg)) ')']);
    end
    byCulture_avg_all(i_cultureID) = nanmean(byCell_avg);
    byCulture_avg_good(i_cultureID) = nanmean(byCell_avg(byCell_good==1));
    disp(['overall (for cells with at least ' num2str(n_cutoff) ' events): ' num2str(byCulture_avg_good(i_cultureID)),10]);

    % -- plotting data -- %
    figure(1); subplot(2,1,1); hold on;
    plot(x_offset + plotdata_good(:,1)+0.15*randn(size(plotdata_good(:,1))),plotdata_good(:,2),'.','Color',plotcolors{i_cultureID},'MarkerSize',5);
    if ~isempty(plotdata_bad)
        h_bad = plot(x_offset + plotdata_bad(:,1)+0.15*randn(size(plotdata_bad(:,1))),plotdata_bad(:,2),'.','Color',(plotcolors{i_cultureID}+[1 1 1])/2,'MarkerSize',5,'LineWidth',1);
        set(get(get(h_bad,'Annotation'),'LegendInformation'),'IconDisplayStyle','off'); % Exclude line from legend
    end
    
    figure(1); subplot(2,1,2); hold on;
    errorbar(x_offset + byCell_cellID(byCell_good==1),byCell_avg(byCell_good==1),byCell_nanstderr(byCell_good==1),'o','Color',plotcolors{i_cultureID},'MarkerSize',4,'LineWidth',0.75);
    if sum(byCell_good) ~= length(byCell_good)
        h_bad = errorbar(x_offset + byCell_cellID(byCell_good==0),byCell_avg(byCell_good==0),byCell_nanstderr(byCell_good==0),'o','Color',(plotcolors{i_cultureID}+[1 1 1])/2,'MarkerSize',4,'LineWidth',0.75);
        set(get(get(h_bad,'Annotation'),'LegendInformation'),'IconDisplayStyle','off'); % Exclude line from legend
    end
    
    x_offset = 1 + x_offset + max(plotdata_all(:,1));
    byCell_ticklabels{end+1} = '';
end

% -- tweek plots -- %
figure(1); subplot(2,1,1)
legend(culturelabels,'Location','NorthEast','Orientation','Vertical')
ylimits = get(gca,'ylim');
axis([0 x_offset 0 ylimits(2)]);
set(gca,'XTick',[0:1:x_offset],'XTickLabel',byCell_ticklabels,'TickDir','out','TickLength',[0.005 0.025]);
ylabel([f_label ' (' f_unit ')'],'FontSize',14);

figure(1); subplot(2,1,2);
%legend(culturelabels,'Location','NorthEast','Orientation','Vertical')
legend(culturelabels,'Location','SouthEast','Orientation','Horizonatal')
ylimits = get(gca,'ylim');
axis([0 x_offset 0 ylimits(2)]);
set(gca,'XTick',[0:1:x_offset],'XTickLabel',byCell_ticklabels,'TickDir','out','TickLength',[0.005 0.025]);
ylabel([f_label ' (' f_unit ')'],'FontSize',14);
xlabel('Cell Number','FontSize',14);

for i_plotmeans = 1:length(cultureID)
    figure(1); subplot(2,1,1)
    plot([0.5 x_offset-0.5],[byCulture_avg_good(i_plotmeans) byCulture_avg_good(i_plotmeans)],'--','Color',plotcolors{i_plotmeans});
    figure(1); subplot(2,1,2);
    plot([0.5 x_offset-0.5],[byCulture_avg_good(i_plotmeans) byCulture_avg_good(i_plotmeans)],'--','Color',plotcolors{i_plotmeans});
end

%% plot by culture %%

% close all;
% f2 = figure(2);
% set(f2,'Position',[100 200 500 400])
% f3 = figure(3);
% set(f3,'Position',[600 200 500 400])
% 

for i_cultureID = 1:length(cultureID)
    load([dirID cultureID{i_cultureID} '_alldata.mat']);
    sampdata = [];
    clear byCulture*
    
    % -- generating sample data -- %
    for i_dat = 1:length(dat)
        if length(dat{i_dat}) >= n_cutoff
            sampdata = [sampdata; randsample(dat{i_dat}(:,f_datcol),n_cutoff)];
        end
    end
    
    bin = 4:1:150;
    byCulture_nhist = histc(sampdata,bin);
    byCulture_cumhist = cumsum(byCulture_nhist);
    byCulture_avg = nanmean(sampdata);
    byCulture_n = length(sampdata);
    byCulture_stderr = nanstd(sampdata)/sqrt(length(sampdata));
    
figure(2); hold on
bar(i_cultureID,byCulture_avg,'FaceColor',plotcolors{i_cultureID},'LineWidth',1)
errorbar(i_cultureID,byCulture_avg,byCulture_stderr,'-k','LineWidth',2);
figure(3); hold on;
plot(1:length(byCulture_cumhist),byCulture_cumhist/max(byCulture_cumhist),'Color',plotcolors{i_cultureID},'LineWidth',3);

% %histograms
% figure(5); hold on
% plot(bin,smooth(byCulture_nhist,5)/max(smooth(byCulture_nhist,5)),'-','Color',plotcolors{i_cultureID},'LineWidth',3);
% axis([0 30 0 1.05]);

% for statistics
if i_cultureID==1
    ctrl_avg = sampdata;
elseif i_cultureID==2
    CNQX_avg = sampdata;
else
    CNQXstim_avg = sampdata;
end

end

figure(2);
set(gca,'Xtick',[1 2 3],'XTickLabel',culturelabels,'TickDir','out','TickLength',[0.005 0.025]);
ylabel([f_label ' (' f_unit ')'],'FontSize',14);
axis square

figure(3);
xlim([0 50]);
set(gca,'TickDir','out','TickLength',[0.005 0.025]);
ylabel('Cumulative Probability','FontSize',14);
xlabel([f_label ' (' f_unit ')'],'FontSize',14);
legend(culturelabels,'Location','SouthEast')
axis square

% for statistics
[h1 p1] = ttest2(ctrl_avg,CNQX_avg)
[h2 p2] = ttest2(CNQX_avg,CNQXstim_avg)
[h3 p3] = ttest2(ctrl_avg,CNQXstim_avg)

% figure(4); hold off;
% plot(ones(size(ctrl_avg))+0.05*randn(size(ctrl_avg)),ctrl_avg,'.','Color',plotcolors{1}); hold on;
% plot(2*ones(size(CNQX_avg))+0.05*randn(size(CNQX_avg)),CNQX_avg,'.','Color',plotcolors{2});
% plot(3*ones(size(CNQXstim_avg))+0.05*randn(size(CNQXstim_avg)),CNQXstim_avg,'.','Color',plotcolors{3});
% axis([0 4 0 40]);
% set(gca,'TickDir','out','TickLength',[0.005 0.025]);



%% save figures %%
figure(1);
export_fig([savename 'ByCell_' num2str(n_cutoff) 'events.png'],'-png','-nocrop','-transparent');
saveas(gcf,[savename 'ByCell_' num2str(n_cutoff) 'events.fig'],'fig')

figure(2);
export_fig([savename 'ByCulture_' num2str(n_cutoff) 'events.png'],'-png','-nocrop','-transparent');
saveas(gcf,[savename 'ByCulture_' num2str(n_cutoff) 'events.fig'],'fig')

if f_cumhist == 1
    figure(3);
    export_fig([savename 'CumHist_' num2str(n_cutoff) 'events.png'],'-png','-nocrop','-transparent');
    saveas(gcf,[savename 'CumHist_' num2str(n_cutoff) 'events.fig'],'fig')
end
