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

%% EXTRACT ORIGINAL DATA %%

plot_title = '';
% load('C:\Users\mingfaifong\Documents\2012-09-13_TTX-CTZ_patch-analysis\ctrl_alldata'); dat_ctrl = dat;
% load('C:\Users\mingfaifong\Documents\2012-09-13_TTX-CTZ_patch-analysis\TTX_alldata'); dat_TTX = dat;
% load('C:\Users\mingfaifong\Documents\2012-09-13_TTX-CTZ_patch-analysis\TTXCTZ_alldata'); dat_TTXCTZ = dat;
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_TTX = dat;
load('C:\Users\mingfaifong\Documents\2012-09-13_CNQX-ChR2_analysis\CNQXChR2_alldata'); dat_TTXCTZ = dat;
% for i = 1:max([size(dat_TTX,1),size(dat_TTXCTZ,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_TTX)
    x_TTX(:,i) = histc(dat_TTX{i}(:,3),edges);
    y_TTX(:,i) = x_TTX(:,i)/size(dat_TTX{i},1);
end

for i = 1:length(dat_TTXCTZ)
    x_TTXCTZ(:,i) = histc(dat_TTXCTZ{i}(:,3),edges);
    y_TTXCTZ(:,i) = x_TTXCTZ(:,i)/size(dat_TTXCTZ{i},1);
end

%% GENERATING SAMPLES DATA%%

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

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

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

maxmini = max([sampdata_ctrl; sampdata_TTX; sampdata_TTXCTZ]);


%% GENERATING DATA FOR CUMULATIVE HISTOGRAM FROM SAMPLE DATA%%
bin = 4:0.5:ceil(maxmini);
cumhist_ctrl = cumsum(histc(sampdata_ctrl,bin));
cumhist_TTX = cumsum(histc(sampdata_TTX,bin));
cumhist_TTXCTZ = cumsum(histc(sampdata_TTXCTZ,bin));


%% GENERATING DATA FOR LINEAR RANK FROM SAMPLE DATA%%
maxsamp = min([length(sampdata_ctrl),length(sampdata_TTX),length(sampdata_TTXCTZ)]);
rankdata_ctrl = randsample(sampdata_ctrl,maxsamp);
rankdata_TTX = randsample(sampdata_TTX,maxsamp);
rankdata_TTXCTZ = randsample(sampdata_TTXCTZ,maxsamp);

rankdata_ctrl = sort(rankdata_ctrl);
rankdata_TTX = sort(rankdata_TTX);
rankdata_TTXCTZ = sort(rankdata_TTXCTZ);

%maxmini = max([rankdata_ctrl; rankdata_TTX; rankdata_TTXCTZ]);

%% 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_TTX,'stacked');
xlim([0 35]); grid on;

subplot(4,1,3);  hold off;
bar(edges+0.25,y_TTXCTZ,'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_TTX,2))/max(smooth(sum(y_TTX,2))),'-r','LineWidth',2);
plot(edges+0.25,smooth(sum(y_TTXCTZ,2))/max(smooth(sum(y_TTXCTZ,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_TTX,'ko');
hold on; plot([0 100],[0 100],'k:');
tmp = regstats(rankdata_TTX,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 TTX mEPSC Amplitude (pA)','FontSize',12);
 
subplot(1,3,2); hold off;
plot(rankdata_ctrl,rankdata_TTXCTZ,'ko');
hold on; plot([0 100],[0 100],'k:');
tmp = regstats(rankdata_TTXCTZ,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 TTX+CTZ mEPSC Amplitude (pA)','FontSize',12);

subplot(1,3,3); hold off;
plot(rankdata_TTX,rankdata_TTXCTZ,'ko');
hold on; plot([0 100],[0 100],'k:');
tmp = regstats(rankdata_TTXCTZ,rankdata_TTX);
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 TTX mEPSC Amplitude (pA)','FontSize',12);
ylabel('Ranked TTX+CTZ 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_TTX/max(cumhist_TTX),'-b','LineWidth',3);
plot(bin,cumhist_TTXCTZ/max(cumhist_TTXCTZ),'-r','LineWidth',3);
legend('control','TTX','TTX+CTZ','Location','SouthEast');
axis square
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_TTX));
errorbar(2,mean(sampdata_TTX),std(sampdata_TTX)/sqrt(length(sampdata_TTX)),'k','LineWidth',2);
bar(3,mean(sampdata_TTXCTZ));
errorbar(3,mean(sampdata_TTXCTZ),std(sampdata_TTXCTZ)/sqrt(length(sampdata_TTXCTZ)),'k','LineWidth',2);
xlim([0.25 3.75]);
set(gca,'XTick',[1 2 3],'XTickLabel',{'control','TTX','TTX+CTZ'},'FontSize',16);
box off;