clear all;

dirID = 'C:\Users\mingfaifong\Documents\2013-01-16_TTX-CTZ-vehicle_patch-analysis\';
type = '30events_01';
ctrl = load([dirID,type,'\ctrl_',type(1:2),'events']);
TTX = load([dirID,type,'\TTX_',type(1:2),'events']);
TTXCTZ = load([dirID,type,'\TTX-CTZ_',type(1:2),'events']);
% ctrl = load([dirID,type,'\ctrl_30events']);
% TTX = load([dirID,type,'\TTX_30events']);
% TTXCTZ = load([dirID,type,'\TTX-CTZ_30events']);

samples.ctrl = reshape(ctrl.sample.mEPSC,1,size(ctrl.sample.mEPSC,1)*size(ctrl.sample.mEPSC,2));
samples.TTX = reshape(TTX.sample.mEPSC,1,size(TTX.sample.mEPSC,1)*size(TTX.sample.mEPSC,2));
samples.TTXCTZ = reshape(TTXCTZ.sample.mEPSC,1,size(TTXCTZ.sample.mEPSC,1)*size(TTXCTZ.sample.mEPSC,2));
maxmini = max([samples.ctrl, samples.TTX, samples.TTXCTZ]);

n_resize = min([length(samples.ctrl),length(samples.TTX),length(samples.TTXCTZ)]);
rankdata.ctrl = sort(randsample(samples.ctrl,n_resize));
rankdata.TTX = sort(randsample(samples.TTX,n_resize));
rankdata.TTXCTZ = sort(randsample(samples.TTXCTZ,n_resize));
stats.ctrl_TTX = regstats(rankdata.TTX,rankdata.ctrl);
stats.TTX_TTXCTZ = regstats(rankdata.TTXCTZ,rankdata.TTX);
stats.ctrl_TTXCTZ = regstats(rankdata.TTXCTZ,rankdata.ctrl);

dt = 1; %pA
bin = 0:dt:ceil(maxmini/dt)*dt;
% f_scale = 1/stats.ctrl_TTX.beta(2);
f_scale = 0.75

x1 = histc(samples.ctrl,bin);
x2 = histc(samples.TTX,bin);

for i = 1:10000
samples.scaledTTX = samples.TTX*f_scale;
samples.scaledTTX = samples.scaledTTX(samples.scaledTTX>=5);
x2_scaled = histc(samples.scaledTTX,bin);

figure(5);
subplot(1,3,1);
hold off;
plot(bin,cumsum(x1)/max(cumsum(x1)),'-k','LineWidth',3);
hold on;
plot(bin,cumsum(x2)/max(cumsum(x2)),'-b','LineWidth',3);
plot(bin,cumsum(x2_scaled)/max(cumsum(x2_scaled)),':b','LineWidth',3);
set(gca,'YTick',[0:0.2:100],'XTick',[0:10:50],'FontSize',14);
xlabel('Amplitude (pA)','FontSize',16)
ylabel('Fraction of Events');
box off;
xlim([0 50]);
ylim([0 1.02]);
axis square;

[h p] = kstest2(x1,x2_scaled)
subplot(1,3,2);
plot(f_scale,p,'ko');
hold on;
axis square

n_resize = min([length(samples.ctrl),length(samples.scaledTTX)]);
rankdata.ctrl_tmp = sort(randsample(samples.ctrl,n_resize));
rankdata.scaledTTX = sort(randsample(samples.scaledTTX,n_resize));
stats.ctrl_scaledTTX = regstats(rankdata.scaledTTX,rankdata.ctrl_tmp);

subplot(1,3,3);
plot(rankdata.ctrl_tmp,rankdata.scaledTTX,'ko');
hold on;
plot([0 maxmini],[0 maxmini],':','LineWidth',1,'Color',[0.5 0.5 0.5]);
hold off
axis square

if round(stats.ctrl_scaledTTX.beta(2)*1000)/1000 == 1
    break;
elseif stats.ctrl_scaledTTX.beta(2) > 1
    f_scale = f_scale*0.9
else
    f_scale = f_scale*1.1
end

pause;
end



%% using middle 90% of distribution

clear all;

dirID = 'C:\Users\mingfaifong\Documents\2013-01-16_TTX-CTZ-vehicle_patch-analysis\';
type = '30events_01';
ctrl = load([dirID,type,'\ctrl_',type(1:2),'events']);
TTX = load([dirID,type,'\TTX_',type(1:2),'events']);
TTXCTZ = load([dirID,type,'\TTX-CTZ_',type(1:2),'events']);
% ctrl = load([dirID,type,'\ctrl_30events']);
% TTX = load([dirID,type,'\TTX_30events']);
% TTXCTZ = load([dirID,type,'\TTX-CTZ_30events']);

samples.ctrl = reshape(ctrl.sample.mEPSC,1,size(ctrl.sample.mEPSC,1)*size(ctrl.sample.mEPSC,2));
samples.TTX = reshape(TTX.sample.mEPSC,1,size(TTX.sample.mEPSC,1)*size(TTX.sample.mEPSC,2));
samples.TTXCTZ = reshape(TTXCTZ.sample.mEPSC,1,size(TTXCTZ.sample.mEPSC,1)*size(TTXCTZ.sample.mEPSC,2));
maxmini = max([samples.ctrl, samples.TTX, samples.TTXCTZ]);

dt = 1; %pA
bin = 0:dt:ceil(maxmini/dt)*dt;

n_resize = min([length(samples.ctrl),length(samples.TTX),length(samples.TTXCTZ)]);
rankdata.ctrl = sort(randsample(samples.ctrl,n_resize));
rankdata.TTX = sort(randsample(samples.TTX,n_resize));
rankdata.TTXCTZ = sort(randsample(samples.TTXCTZ,n_resize));

range90 = [1:n_resize];
%range90 = [round(n_resize*0.4):round(n_resize*0.6)];
stats.ctrl_TTX = regstats(rankdata.TTX(range90),rankdata.ctrl(range90)); stats.ctrl_TTX.beta
stats.TTX_TTXCTZ = regstats(rankdata.TTXCTZ(range90),rankdata.TTX(range90));
stats.ctrl_TTXCTZ = regstats(rankdata.TTXCTZ(range90),rankdata.ctrl(range90));

x1 = histc(samples.ctrl,bin);
x2 = histc(samples.TTX,bin);
samples.scaledTTX = (samples.TTX-stats.ctrl_TTX.beta(1))/stats.ctrl_TTX.beta(2);
%samples.scaledTTX = samples.scaledTTX(samples.scaledTTX>=5);
x2_scaled = histc(samples.scaledTTX,bin);
x3 = histc(samples.TTXCTZ,bin);
samples.scaledTTXCTZ = (samples.TTXCTZ-stats.ctrl_TTXCTZ.beta(1))/stats.ctrl_TTXCTZ.beta(2);
%samples.scaledTTXCTZ = samples.scaledTTXCTZ(samples.scaledTTXCTZ>=5);
x3_scaled = histc(samples.scaledTTXCTZ,bin);

figure(6);
hold off;
plot(bin,cumsum(x1)/max(cumsum(x1)),'-k','LineWidth',3);
hold on;
plot(bin,cumsum(x2)/max(cumsum(x2)),'-b','LineWidth',3);
plot(bin,cumsum(x2_scaled)/max(cumsum(x2_scaled)),':b','LineWidth',3);
plot(bin,cumsum(x3)/max(cumsum(x3)),'-r','LineWidth',3);
plot(bin,cumsum(x3_scaled)/max(cumsum(x3_scaled)),':r','LineWidth',3);
set(gca,'YTick',[0:0.2:100],'XTick',[0:10:50],'FontSize',14);
xlabel('Amplitude (pA)','FontSize',16)
ylabel('Fraction of Events');
box off;
xlim([0 50]);
ylim([0 1.02]);
axis square;

[h1 p1] = kstest2(x1,x2_scaled)
[h1a p1a] = ranksum(x1,x2_scaled)

[h2 p2] = kstest2(x2_scaled,x3_scaled)
[h2a p2a] = ranksum(x2_scaled,x3_scaled)

[h3 p3] = kstest2(x1,x3_scaled)
[h3a p3a] = ranksum(x1,x3_scaled)

% 
% x4 = histc(samples.TTXCTZ,bin);
% samples.scaledTTXCTZ = (samples.TTXCTZ-stats.TTX_TTXCTZ.beta(1))/stats.TTX_TTXCTZ.beta(2);
% %samples.scaledTTXCTZ = samples.scaledTTXCTZ(samples.scaledTTXCTZ>=5);
% x4_scaled = histc(samples.scaledTTXCTZ,bin);
% [h4 p4] = kstest2(x2,x4_scaled)
% [h4a p4a] = ranksum(x2,x4_scaled)
% 
% plot(bin,cumsum(x4_scaled)/max(cumsum(x4_scaled)),':r','LineWidth',3);