clear all;
f_datcol = 3; %amplitudes
n_cutoff = 19;
isctrl = 0;

%% EXTRACT ORIGINAL DATA %%

plot_title = '';
load('C:\Users\mingfaifong\Documents\2012-09-13_CNQX-ChR2_analysis\ctrl_alldata'); dat_ctrl = dat;
load('C:\Users\mingfaifong\Documents\2012-09-13_CNQX-ChR2_analysis\CNQX_alldata'); dat_CNQX = dat;
load('C:\Users\mingfaifong\Documents\2012-09-13_CNQX-ChR2_analysis\CNQXChR2_alldata'); dat_CNQXChR2 = dat;
% for i = 1:max([size(dat_CNQX,1),size(dat_CNQXChR2,1)])
%     dat_ctrl{i} = nan(100,19);
% end


%% GENERATING DATA FOR AMPLITUDE HISTOGRAM FROM ALL DATA%%
edges = 0:1:300;
for i = 1:length(dat_ctrl)
    x_ctrl(:,i) = histc(dat_ctrl{i}(:,3),edges);
    y_ctrl(:,i) = x_ctrl(:,i)/size(dat_ctrl{i},1);
end

for i = 1:length(dat_CNQX)
    x_CNQX(:,i) = histc(dat_CNQX{i}(:,3),edges);
    y_CNQX(:,i) = x_CNQX(:,i)/size(dat_CNQX{i},1);
end

for i = 1:length(dat_CNQXChR2)
    x_CNQXChR2(:,i) = histc(dat_CNQXChR2{i}(:,3),edges);
    y_CNQXChR2(:,i) = x_CNQXChR2(:,i)/size(dat_CNQXChR2{i},1);
end

%% GENERATING SAMPLES DATA%%

sampdata_ctrl = [];
for i = 1:size(dat_ctrl,1)
%    if length(dat_ctrl{i}) >= n_cutoff
    if size(dat_ctrl{i},1) >= n_cutoff
        sampdata_ctrl = [sampdata_ctrl; randsample(dat_ctrl{i}(:,f_datcol),n_cutoff)];
    end
end

sampdata_CNQX = [];
for i = 1:size(dat_CNQX,1)
    %if length(dat_CNQX{i}) >= n_cutoff
    if size(dat_CNQX{i},1) >= n_cutoff
        sampdata_CNQX = [sampdata_CNQX; randsample(dat_CNQX{i}(:,f_datcol),n_cutoff)];
    end
end

sampdata_CNQXChR2 = [];
for i = 1:size(dat_CNQXChR2,1)
    %if length(dat_CNQXChR2{i}) >= n_cutoff
    if size(dat_CNQXChR2{i},1) >= n_cutoff
        sampdata_CNQXChR2 = [sampdata_CNQXChR2; randsample(dat_CNQXChR2{i}(:,f_datcol),n_cutoff)];
    end
end

maxmini = max([sampdata_ctrl; sampdata_CNQX; sampdata_CNQXChR2]);


%% GENERATING DATA FOR CUMULATIVE HISTOGRAM FROM SAMPLE DATA%%
bin = 4:0.5:ceil(maxmini);
cumhist_ctrl = cumsum(histc(sampdata_ctrl,bin));
cumhist_CNQX = cumsum(histc(sampdata_CNQX,bin));
cumhist_CNQXChR2 = cumsum(histc(sampdata_CNQXChR2,bin));


%% GENERATING DATA FOR LINEAR RANK FROM SAMPLE DATA%%
maxsamp = min([length(sampdata_ctrl),length(sampdata_CNQX),length(sampdata_CNQXChR2)]);
rankdata_ctrl = randsample(sampdata_ctrl,maxsamp);
rankdata_CNQX = randsample(sampdata_CNQX,maxsamp);
rankdata_CNQXChR2 = randsample(sampdata_CNQXChR2,maxsamp);

rankdata_ctrl = sort(rankdata_ctrl);
rankdata_CNQX = sort(rankdata_CNQX);
rankdata_CNQXChR2 = sort(rankdata_CNQXChR2);

%maxmini = max([rankdata_ctrl; rankdata_CNQX; rankdata_CNQXChR2]);

%% PLOTTING %%

figure(1);
subplot(4,1,1);  hold off;
bar(edges+0.25,y_ctrl,'stacked');
xlim([0 35]); grid on;

subplot(4,1,2);  hold off;
bar(edges+0.25,y_CNQX,'stacked');
xlim([0 35]); grid on;

subplot(4,1,3);  hold off;
bar(edges+0.25,y_CNQXChR2,'stacked');
xlim([0 35]); grid on;

subplot(4,1,4); hold off;
plot(edges+0.25,smooth(sum(y_ctrl,2))/max(smooth(sum(y_ctrl,2))),'-k','LineWidth',2);
hold on;
plot(edges+0.25,smooth(sum(y_CNQX,2))/max(smooth(sum(y_CNQX,2))),'-r','LineWidth',2);
plot(edges+0.25,smooth(sum(y_CNQXChR2,2))/max(smooth(sum(y_CNQXChR2,2))),'-b','LineWidth',2);
xlim([0 35]); grid on;
ylim([0 1.1]);


figure(2);
subplot(1,3,1); hold off;
plot(rankdata_ctrl,rankdata_CNQX,'ko');
hold on; plot([0 100],[0 100],'k:');
tmp = regstats(rankdata_CNQX,rankdata_ctrl);
plot([0 100],tmp.beta(1)+[0 100]*tmp.beta(2),'-k','LineWidth',2);
text(5,95,['r^2 = ' num2str(tmp.rsquare)],'FontSize',12);
axis([0 100 0 100]);
axis square; box off;
xlabel('Ranked Ctrl mEPSC Amplitude (pA)','FontSize',12);
ylabel('Ranked CNQX mEPSC Amplitude (pA)','FontSize',12);
 
subplot(1,3,2); hold off;
plot(rankdata_ctrl,rankdata_CNQXChR2,'ko');
hold on; plot([0 100],[0 100],'k:');
tmp = regstats(rankdata_CNQXChR2,rankdata_ctrl);
plot([0 100],tmp.beta(1)+[0 100]*tmp.beta(2),'-k','LineWidth',2);
text(5,95,['r^2 = ' num2str(tmp.rsquare)],'FontSize',12);
axis([0 100 0 100]);
axis square; box off;
xlabel('Ranked Ctrl mEPSC Amplitude (pA)','FontSize',12);
ylabel('Ranked CNQX+stim mEPSC Amplitude (pA)','FontSize',12);

subplot(1,3,3); hold off;
plot(rankdata_CNQX,rankdata_CNQXChR2,'ko');
hold on; plot([0 100],[0 100],'k:');
tmp = regstats(rankdata_CNQXChR2,rankdata_CNQX);
plot([0 100],tmp.beta(1)+[0 100]*tmp.beta(2),'-k','LineWidth',2);
text(5,95,['r^2 = ' num2str(tmp.rsquare)],'FontSize',12);
axis([0 100 0 100]);
axis square; box off; 
xlabel('Ranked CNQX mEPSC Amplitude (pA)','FontSize',12);
ylabel('Ranked CNQX+stim mEPSC Amplitude (pA)','FontSize',12);

figure(3); hold off;
bin = 4:0.5:ceil(maxmini);
plot(bin,cumhist_ctrl/max(cumhist_ctrl),'-k','LineWidth',3);
hold on;
plot(bin,cumhist_CNQX/max(cumhist_CNQX),'-b','LineWidth',3);
plot(bin,cumhist_CNQXChR2/max(cumhist_CNQXChR2),'-r','LineWidth',3);
legend('control','CNQX','CNQX+stim','Location','SouthEast');
axis square; box off;
axis([0 75 0 1.05]);
xlabel('mEPSC Amplitude (pA)','FontSize',16); ylabel('Fraction of Events','FontSize',16);


figure(4); hold off;
bar(1,mean(sampdata_ctrl));
hold on;
errorbar(1,mean(sampdata_ctrl),std(sampdata_ctrl)/sqrt(length(sampdata_ctrl)),'k','LineWidth',2);
bar(2,mean(sampdata_CNQX));
errorbar(2,mean(sampdata_CNQX),std(sampdata_CNQX)/sqrt(length(sampdata_CNQX)),'k','LineWidth',2);
bar(3,mean(sampdata_CNQXChR2));
errorbar(3,mean(sampdata_CNQXChR2),std(sampdata_CNQXChR2)/sqrt(length(sampdata_CNQXChR2)),'k','LineWidth',2);
xlim([0.25 3.75]);
set(gca,'XTick',[1 2 3],'XTickLabel',{'control','CNQX','CNQX+stim'},'FontSize',16);
box off;
axis square;