% % Jon Newman
% 
% %% file id
% fids = {
%     'J:/newmonster/2012_02_09_Lollipop_Long-Term-optoBB/cl-long-bb_01_2012-02-09_chr2_lollipop_10803';      % 0.75 Hz
%     'J:/newmonster/2012_02_10_Lollipop_Long-Term-optoBB/cl-long-bb_01_2012-02-10_chr2_lollipop_10803_1Hz';  % 1 Hz
%     'J:/newmonster/2012_02_27_Lollipop_Long-Term-optoBB/2012_02_27_Lollipop_Long-Term-optoBB_00';           % 2 Hz
%     'J:/newmonster/2012_02_11_Lollipop_Long-Term-optoBB/cl-long-bb_01_2012-02-11_chr2_lollipop_10803';      % 3 Hz
%     'J:/newmonster/2012_02_16_Lollipop_Long-Term-optoBB/cl-long-bb_01_2012-02-16_chr2_lollipop_10803';      % 4 Hz
%     'J:/newmonster/2012_02_14_Lollipop_Long-Term-optoBB/cl-long-bb_01_2012-02-14_chr2_lollipop_10803';      % 5 Hz
%     'J:/newmonster/2012_02_12_Lollipop_Long-Term-optoBB/cl-long-bb_01_2012-02-12_chr2_lollipop_10803';      % 6 Hz
%     };
% 
% 
% %% Parameters
% % time bounds for the plots
% tbound = [0 24*3600-60]; %3 hrs pre, 12 hrs clamp, 9 hrs post.
% sort_bound = [0 3*3600]; % The time to sort unit firing rate over
% 
DIV = [9 10 27 11 16 15 13] - 9 + 70;
HRS = [20 20 20 23 20 15 16];

% 
% % Upper bound for UNFR plot
% maxFR = 10;
% 
% % regions of interests for zoomed rasters, in hours
% r1 = 1.5;
% r2 = 9.0;
% r3 = 16;
% r4 = 24;
% roiwidth = 60; % region of interest width in seconds
% s_T = 10;
% tdev = 3600; % time (in seconds) scaling factor
% 
% % drug parameters
% isdrug = 0;
% tdrug = 3600*8;
% twash = 3600*10;
% 
% dT = 300; % binsize seconds, for all time seris analysis
% stimHz = 10;
% stimWidthSec = 0.005;
% 
% col =  [[100 100 100];
%     [220 000 000];
%     [255 102 000];
%     [255 204 000];
%     [000 170 000];
%     [000 170 212];
%     [085 000 212]]/255;
% 
% 
% %% Figures
% close all
% 
% %
% % Axis layout:
% %
% % [A1a            ]
% % [A1b            ]
% % [A1c            ]
% % [A1d            ]
% % [A1e            ]
% % [A1f            ]
% %
% % |---------------|
% % |A2             |
% % |---------------|
% %
% % |---------------|
% % |A3             |
% % |---------------|
% 
% % Heights and widths
% h1x = 0.65;
% h2  = 1.8;
% h3  = 1.8;
% w   = 11.5;
% w2 = w *( 18/24);
% s   = 0.35; % verticle spacing
% s2  = 0.15;
% s3 = 0.75; % secondary horizontal spacing
% w3 = w *( 6/24) - s3;
% 
% 
% numsp = length(fids); % number of subplots
% percent = 0;
% 
% targets = [];
% uFRpre = [];
% uFRpost = [];
% uFRclamp = [];
% uStimFreq = [];
% medStimFreq = [];
% minStimFreq = [];
% maxStimFreq = [];
% 
% 
% %%
% for f = 1:length(fids)
%     
%     % Assing file
%     fid = fids{f};
%     
%     % Load simoc file
%     sim = loadsimoc([fid '.simoc']);
%     
%     % Load SS
%     load([fid '.SS'],'-mat');
%     
%     % Recalculate ASDR
%     SS.ASDR(dT,[],[],[],[],0);
%     
%     % target values
%     iTargOn = find(sim.target > 0,1,'first');
%     iTargOff = find(sim.target > 0,1,'last');
%     targetVal = sim.target(iTargOn);
%     if targetVal == 1.5
%         targetVal = 1.0;
%     end
%     targets = [targets, targetVal];
%     tTargOn = sim.time(iTargOn);
%     tTargOff = sim.time(iTargOff);
%     
%     %% Mean firing rate before,during, and after clamp
%     uFRpre = [uFRpre numel(SS.time(SS.unit~=0 & SS.time < tTargOn))/numel(unique(SS.unit))/tTargOn];
%     uFRclamp = [uFRclamp numel(SS.time(SS.unit~=0 & SS.time > tTargOn & SS.time <= tTargOff))/numel(unique(SS.unit))/(tTargOff - tTargOn)];
%     uFRpost = [uFRpost numel(SS.time(SS.unit~=0 & SS.time > tTargOff & SS.time <= tbound(2)))/numel(unique(SS.unit))/(tbound(2) - tTargOff)];
%     
%     %% bang-bang input signal (duty cycle)
% %     maxduty = stimWidthSec*stimHz;
%     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)
%     
%     uStimFreq = [uStimFreq mean(stimFreq(stimFreq~=0))];
%     medStimFreq = [medStimFreq median(stimFreq(stimFreq~=0))];
%     minStimFreq = [minStimFreq min(stimFreq(stimFreq~=0))];
%     maxStimFreq = [maxStimFreq max(stimFreq(stimFreq~=0))];
%     
% 
%     
% end



%%

DIE = DIV - min(DIV) + 1 + (HRS- min(HRS))/24;

close all
figure('Units','Centimeter','Position',[5 5 12 12])
subplot(221)

hold on
[~, ind] = sort(DIE);
plot(DIE(ind), uFRpre(ind) ,'ko-','MarkerSize',4)
xlabel('Days into experiment','FontSize',8,'FontName','Arial')
ylabel('FR_{pre} (Hz/unit)','FontSize',8,'FontName','Arial')
axis square
set(gca,'TickDir','out','TickLength',[0.025 0.025],'FontSize',8,'FontName','Arial')

% Linear regression
x = DIE; y = uFRpre;
p = polyfit(x,y,1);
yfit = polyval(p,x);
yresid = y - yfit;
SSresid = sum(yresid.^2);
SStotal = (length(y)-1) * var(y);
rsq = 1 - SSresid/SStotal;
title(['R^{2} = ' num2str(rsq)])
% plot(x,yfit)

subplot(223)

hold on
plot(uFRpre, uFRpost ,'ko','MarkerSize',4)
xlabel('FR_{pre} (Hz/unit)','FontSize',8,'FontName','Arial')
ylabel('FR_{post} (Hz/unit)','FontSize',8,'FontName','Arial')
axis square
set(gca,'TickDir','out','TickLength',[0.025 0.025],'Xlim',[0 3],'Ylim',[0 3],'FontSize',8,'FontName','Arial')

% Linear regression
x = uFRpre; y = uFRpost;
p = polyfit(x,y,1);
yfit = polyval(p,x);
yresid = y - yfit;
SSresid = sum(yresid.^2);
SStotal = (length(y)-1) * var(y);
rsq = 1 - SSresid/SStotal;
title(['R^{2} = ' num2str(rsq)])
% plot(x,yfit,'k-')
plot([0 3],[0 3],'k:')

subplot(224)

hold on
plot(uFRpre, uStimFreq ,'ko','MarkerSize',4)
axis square
xlabel('FR_{pre} (Hz/unit)','FontSize',8,'FontName','Arial')
ylabel('Avg. stim. freq. (Hz)','FontSize',8,'FontName','Arial')
set(gca,'TickDir','out','TickLength',[0.025 0.025],'Xlim',[0 3],'FontSize',8,'FontName','Arial')

% Linear regression
x = uFRpre; y = uStimFreq;
p = polyfit(x,y,1);
yfit = polyval(p,x);
yresid = y - yfit;
SSresid = sum(yresid.^2);
SStotal = (length(y)-1) * var(y);
rsq = 1 - SSresid/SStotal;
title(['R^{2} = ' num2str(rsq)])
% plot(x,yfit)

subplot(222)

hold on
plot(uFRpre, uStimFreq(ind) ,'ko-','MarkerSize',4)
xlabel('Days into experiment','FontSize',8,'FontName','Arial')
ylabel('Avg. stim. freq. (Hz)','FontSize',8,'FontName','Arial')
axis square
set(gca,'TickDir','out','TickLength',[0.025 0.025],'FontSize',8,'FontName','Arial')

% Linear regression
x = uFRpre; y = uStimFreq;
p = polyfit(x,y,1);
yfit = polyval(p,x);
yresid = y - yfit;
SSresid = sum(yresid.^2);
SStotal = (length(y)-1) * var(y);
rsq = 1 - SSresid/SStotal;
title(['R^{2} = ' num2str(rsq)])
% plot(x,yfit)



%% Export figures

export_fig 'stats-summary_cl-long-bb_lollipop_10803_summary_01' -pdf -transparent
saveas(gcf,'stats-summary_cl-long-bb_lollipop_10803_summary_01.fig')
