% % Jon Newman Combine disparate data sets to create large summary figure for
% % optoclamp paper. Because the data are in different formats and have
% % different experimental nauances, this script is going to be a complete,
% % one-off mess.
% 
% Start by making axes
close all
figure('units','centimeter','position',[5 5 20 20]);

% Dimensions (cm)
xo = 2;
yo = 2;
w1 = 12.025;
h1 = 0.8125;
s1 = 0.25;
s2 = 0.75;
s3 = 0.75;
s4 = 0.5;
w2 = 3.5;
h3 = 1.75;
w3 = 1.5;
h2 = 1.5;

% Plot limits (hrs)
tgap = 0/60; 
tpre = 3;
tdrug1 = 5/60; % 
tclamp = 24;
tdrug2 = 5/60;
tpost = 6;
ttotal = (tpre + tdrug1 + tclamp + tdrug2 + tpost + 2*tgap);
% nchunk = 10;
% tchunk = 0:ttotal/(nchunk+1):ttotal;
tcorr = 5/60;


% Make Axes
a = zeros(1,6);
a(4) = axes('units','centimeter','position',[xo yo w1 h1]);
a(3) = axes('units','centimeter','position',[xo yo+h1+s1 w1 h1]);
a(2) = axes('units','centimeter','position',[xo yo+2*(h1+s1) w1 h1]);
a(1) = axes('units','centimeter','position',[xo yo+3*(h1+s1) w1 h1]);

a(6) = axes('units','centimeter','position',[xo+s2+w1 yo w2 h2]);
a(5) = axes('units','centimeter','position',[xo+s2+w1 yo+h2+s3 w2 h3]);

 
%% Experiment 1 (20 uM CNQX) : 2012_03_05_Lollipop-10803

% Storage
sp.time = [];
sp.channel = [];
sp.unit = [];

% Load the relavent spike files
x = load('J:\newmonster\2012_03_05_Lollipop-10803_Long-Term-optoBB_CNQX/2012_03_05_Lollipop-10803_Long-Term-optoBB_CNQX_00_pre-drug.mat');

b = x.s.time > max(x.s.time) - 3600*(tpre) & x.s.unit > 0;
fr_pre = sum(b)/(3600*tpre);
sp.time = [sp.time;  x.s.time(b) - (max(x.s.time) - 3600*(tpre))];
sp.channel = [sp.channel;  x.s.channel(b)];
sp.unit = [sp.unit;  x.s.unit(b)];
clear x

sim = loadsimoc('J:\newmonster\2012_03_05_Lollipop-10803_Long-Term-optoBB_CNQX/2012_03_05_Lollipop-10803_Long-Term-optoBB_CNQX_00_20uM-CNQX.simoc');
t_start = min(sim.time(sim.target > 0));
t_end = max(sim.time(sim.target > 0));
tclamp = t_end-t_start;

x = load('J:\newmonster\2012_03_05_Lollipop-10803_Long-Term-optoBB_CNQX/2012_03_05_Lollipop-10803_Long-Term-optoBB_CNQX_00_20uM-CNQX.mat');

b = x.s.time > t_start - tdrug1*3600 &  x.s.time < t_end + 3600*tdrug2 & x.s.unit > 0;
fr_drug1 = sum(x.s.time(b) < t_start)/(tdrug1*3600);
fr_clamp = sum(x.s.time(b) > t_start & x.s.time(b) < t_end)/tclamp;
fr_drug2 = sum(x.s.time(b) > t_end)/(3600*tdrug2);
sp.time = [sp.time;  x.s.time(b) - (t_start - tdrug1*3600) + 3600*tpre + 3600*tgap];
sp.channel = [sp.channel;  x.s.channel(b)];
sp.unit = [sp.unit;  x.s.unit(b)];
clear x
 
x = load('J:\newmonster\2012_03_05_Lollipop-10803_Long-Term-optoBB_CNQX/2012_03_05_Lollipop-10803_Long-Term-optoBB_CNQX_00_wash4-5.mat');

b = x.s.time < 3600*(tpost) & x.s.unit > 0;
fr_post = sum(b)/(3600*tpost);
sp.time = [sp.time;  x.s.time(b) + 3600*tpre + tclamp + 3600*tdrug1 + 3600*tdrug2 + 3600*2*tgap];
sp.channel = [sp.channel;  x.s.channel(b)];
sp.unit = [sp.unit;  x.s.unit(b)];
clear x

ss = SqueakySpk([],25000,1,sp);
clear sp

%%
dt = 1;
ss.ASDR(dt,[],[],[],[],0);


%%
axes(a(2))
hold on
stairs(ss.asdr.bin/3600 - tpre,ss.asdr.asdr/max(ss.unit),'k-')
set(a(2),'box','off','XTick',[],'TickDir','out','Xlim',[-tpre ttotal-tpre],'Ylim',[0 25])

%%
dt = 300;
ss.ASDR(dt,[],[],[],[],0);

%%
axes(a(1))
ss.PlotUSDR(7,1,[0 3600*tpre],0)
set(a(1),'box','off','XTick',[],'TickDir','out')
xlabel('')

%%
axes(a(3))
stairs(ss.asdr.bin/3600 - tpre,ss.asdr.asdr/max(ss.unit),'k-')
hold on
targ = unique(sim.target);
targ = targ(targ>0);
plot([tpre + tgap + tdrug1, tpre + tgap + tclamp/3600] - tpre,[targ, targ],'r-')
set(a(3),'box','off','XTick',[],'TickDir','out','Xlim',[-tpre ttotal-tpre])

%%
axes(a(4))
    tduty = sim.time(sim.output(1,:) > 0);
    c = histc(tduty, ss.asdr.bin);
    stimFreq = stimHz*c/dt/100; % convert to stim freq (update rate is 100 Hz, so need to account for that)

stairs(ss.asdr.bin/3600,stimFreq,'k-')
set(a(4),'box','off','TickDir','out','Xlim',[-tpre ttotal-tpre])

%%
axes(a(5))
hold on
plot([1 2 3 4 5],[fr_pre fr_drug1 fr_clamp fr_drug2 fr_post]/fr_pre,'ko-')
plot([0.5 5.5],[1 1],'k:')
set(a(5),'box','off','TickDir','out')

%%
axes(a(6))
hold on
mpre = tpre/2;
generate_xcorr(ss,3600*[mpre-tcorr, mpre+tcorr])


mdrug1 = tpre + tgap + tdrug1/2;
generate_xcorr(ss,3600*[mdrug1 - tcorr/2,  mdrug1 + tcorr/2])


mclamp = tpre + tgap + tdrug1 + tclamp/2/3600;
generate_xcorr(ss,3600*[mclamp - tcorr,  mclamp + tcorr])


mdrug2 = tpre + tgap + tdrug1 + tclamp/3600 + tdrug2/2;
generate_xcorr(ss,3600*[mdrug2  - tcorr,  mdrug2  + tcorr])

mpost = tpre + tgap + tdrug1 + tclamp/3600 + tdrug2 + tpost/2;
generate_xcorr(ss,3600*[mclamp - tcorr,  mpost + tcorr])

%
hand = get(a(6), 'Children');
set(hand(5),'Color',[0 100 255]/255)
set(hand(4),'Color',[220 0 0]/255)
set(hand(3),'Color',[0 170 0]/255)
set(hand(2),'Color',[155 0 212]/255)
set(hand(1),'Color',[255 204 0]/255)

set(a(6),'box','off','TickDir','out','Ylim',[-0.5 8])

%% Raster plots of various size
wp1 = 5*(w1+w2+s2-3*s1)/16;
ra(1) = axes('units','centimeter','position',[xo yo+4*h1 + 4*s1 wp1 h1/1.75]);
ss.RasterPlot(3600*[mpre-2/3600, mpre+3/3600]) % 10 seconds
set(ra(1),'box','on','TickDir','out','XTick',[],'YTick',[])
xlabel ''
ylabel ''
title ''

axes(a(2))
plot([mpre-2/3600, mpre+3/3600]-tpre,[20 20],'r*')

%
wp2 = 2*(w1+w2+s2-3*s1)/16;
ra(2) = axes('units','centimeter','position',[xo+wp1+s1 yo+4*h1 + 4*s1 wp2 h1/1.75]);
ss.RasterPlot(3600*[mdrug1-1/3600,mdrug1+1/3600])
set(ra(2),'box','on','TickDir','out','XTick',[],'YTick',[])
xlabel ''
ylabel ''
title ''

axes(a(2))
plot([mdrug1-1/3600,mdrug1+1/3600]-tpre,[20 20],'r*')

%
wp3 = 4*(w1+w2+s2-3*s1)/16;
ra(3) = axes('units','centimeter','position',[xo+wp1+wp2+2*s1 yo+4*h1 + 4*s1 wp3 h1/1.75]);
ss.RasterPlot(3600*[mclamp+15/3600,mclamp+19/3600])
set(ra(3),'box','on','TickDir','out','XTick',[],'YTick',[])
xlabel ''
ylabel ''
title ''

axes(a(2))
plot([mclamp+15/3600,mclamp+19/3600]-tpre,[20 20],'r*')


%
wp4 = 5*(w1+w2+s2-3*s1)/16;
ra(4) = axes('units','centimeter','position',[xo+wp1+wp2+wp3+3*s1 yo+4*h1 + 4*s1 wp4 h1/1.75]);
ss.RasterPlot(3600*[mpost+15/3600,mpost+20/3600])
set(ra(4),'box','on','TickDir','out','XTick',[],'YTick',[])
xlabel ''
ylabel ''
title ''

axes(a(2))
plot([mpost+15/3600,mpost+20/3600]-tpre,[20 20],'r*')



%% Experiment 2 (40 uM CNQX) : 2013_02_26_Yam-22565
% Load the relavent spike files

% Storage
sp.time = [];
sp.channel = [];
sp.unit = [];

s = loadspk('Z:\jon\2013_02_26_Yam-22565_24-hr-cnqx-clamp\2013_02_26_Yam-22565_no-drug_01.spk');

b = s.time > max(s.time) - 3600*(tpre) & s.unit > 0;
fr_pre = sum(b)/(3600*tpre);
sp.time = [sp.time;  s.time(b) - (max(s.time) - 3600*(tpre))];
sp.channel = [sp.channel;  s.channel(b)];
sp.unit = [sp.unit;  s.unit(b)];
clear x

bang = loadbang('Z:\jon\2013_02_26_Yam-22565_24-hr-cnqx-clamp\2013_02_26_Yam-22565_40-um-cnqx_01.bang');


t_start = min(bang.time(bang.clon > 0));
t_end = max(bang.time(bang.clon > 0));
tclamp = t_end-t_start;

s = loadspk('Z:\jon\2013_02_26_Yam-22565_24-hr-cnqx-clamp\2013_02_26_Yam-22565_40-um-cnqx_01.spk');

b = s.time > t_start - tdrug1*3600 &  s.time < t_end + 3600*tdrug2 & s.unit > 0;
fr_drug1 = sum(s.time(b) < t_start)/(tdrug1*3600);
fr_clamp = sum(s.time(b) > t_start & s.time(b) < t_end)/tclamp;
fr_drug2 = sum(s.time(b) > t_end)/(3600*tdrug2);
sp.time = [sp.time;  s.time(b) - (t_start - tdrug1*3600) + 3600*tpre + 3600*tgap];
sp.channel = [sp.channel;  s.channel(b)];
sp.unit = [sp.unit;  s.unit(b)];
clear x
 
s = loadspk('Z:\jon\2013_02_26_Yam-22565_24-hr-cnqx-clamp\2013_02_26_Yam-22565_wash7-8_01.spk');

b = s.time < 3600*(tpost) & s.unit > 0;
fr_post = sum(b)/(3600*tpost);
sp.time = [sp.time;  s.time(b) + 3600*tpre + tclamp + 3600*tdrug1 + 3600*tdrug2 + 3600*2*tgap];
sp.channel = [sp.channel;  s.channel(b)];
sp.unit = [sp.unit;  s.unit(b)];
clear x

ss = SqueakySpk([],25000,1,sp);
clear sp

%%
axes(a(5))
hold on
plot([1 2 3 4 5],[fr_pre fr_drug1 fr_clamp fr_drug2 fr_post]/fr_pre,'ro-')


%%
saveas(gcf,'24_hour_cnqx_decoupling_summary_06');
export_fig -pdf -transparent 24_hour_cnqx_decoupling_summary_06.pdf

