% Jon Newman

%% file id
fids = { 
    'J:/newmonster/2012_02_09_Lollipop_Long-Term-optoBB/cl-long-bb_01_2012-02-09_chr2_lollipop_10803';
    'J:/newmonster/2012_02_10_Lollipop_Long-Term-optoBB/cl-long-bb_01_2012-02-10_chr2_lollipop_10803';
    'J:/newmonster/2012_02_11_Lollipop_Long-Term-optoBB/cl-long-bb_01_2012-02-11_chr2_lollipop_10803';
    'J:/newmonster/2012_02_12_Lollipop_Long-Term-optoBB/cl-long-bb_01_2012-02-12_chr2_lollipop_10803';
    'J:/newmonster/2012_02_14_Lollipop_Long-Term-optoBB/cl-long-bb_01_2012-02-14_chr2_lollipop_10803';
    'J:/newmonster/2012_02_16_Lollipop_Long-Term-optoBB/cl-long-bb_01_2012-02-16_chr2_lollipop_10803';
    'J:/newmonster/2012_02_27_Lollipop_Long-Term-optoBB/2012_02_27_Lollipop_Long-Term-optoBB_00';
    };


%% Parameters
% time bounds for the plots
tbound = [0 24*3600]; %3 hrs pre, 12 hrs clamp, 9 hrs post.
sort_bound = [0 3*3600]; % The time to sort unit firing rate over
% 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;

dT1 = 0.10; % short time scale
dT2 = 300; % binsize seconds, for all time seris analysis
stimHz = 10;
stimWidthSec = 0.01;
maxduty = stimWidthSec*stimHz;

%%
col = [[220 000 000];
[255 102 000];
[255 204 000];
[000 170 000];
[000 170 212];
[000 100 255];
[085 000 212]]/255;

% Target UNDR's
targets = [0.75 1.5 2 3 4 5 6];

% Axis positions
close all
f1 = figure('units','centimeter','position',[0 0 18.5 20.25]);
pp1a = [];
pp1b = [];
yo = 10.25;
s = 0.25;
lo = (10-(length(targets)-1)*0.25);
th = lo/length(targets);
th1 = 1*th;
th2 = 0*th;
for i = 1:length(targets)
    pp1a = [pp1a; [1 yo 16.5 th1]];
    pp1b = [pp1b; [1 yo+th1 16.5 th2]];
    yo = yo+th+s;
end

pp2 = [1 5.5 16.5 4];
aclamp = axes('units','centimeters','position',pp2);
pp3 = [1 0.75 16.5 4];
aduty = axes('units','centimeters','position',pp3);


for f = 1:length(fids)
    
    % Assing file
    fid = fids{f};

    % Load simoc file
    try
        sim = loadsimoc([fid '.simoc']);
        procsim = 1;
    catch
        procsim = 0;
        bang = loadbang([fid '.bang']);
    end
        
    % Load SS
    load([fid '.SS'],'-mat');
    
    % Recalculate ASDR
    SS.ASDR(dT2,[],[],[],[],0);
       
    % target values
    if procsim
        
        % targets
        iTargOn = find(sim.target > 0,1,'first');
        iTargOff = find(sim.target > 0,1,'last');
        targetVal = sim.target(iTargOn);
        tTargOn = sim.time(iTargOn);
        tTargOff = sim.time(iTargOff);
        
        % Duty Cycle
        tduty = sim.time(sim.output(1,:) ~= 0);
        c = hist(tduty, SS.asdr.bin);
        duty = maxduty*(c/(100*dT2));
    else
        iTargOn = find(bang.clon > 0,1,'first');
        iTargOff = find(bang.clon  > 0,1,'last');
        targetVal = bang.targ (iTargOn);
        tTargOn = bang.time (iTargOn);
        tTargOff = bang.time(iTargOff);
        
        % Duty Cycle
        tduty = bang.time(bang.out1 ~= 0);
        c = hist(tduty, SS.asdr.bin);
        duty = maxduty*(c/(100*dT2));
    end
    

    %% bang-bang input signal (duty cycle)
    
    %% USDR
    ind = find(targets == targetVal);
    indr = (length(targets) - ind) + 1;
    axes('units','centimeters','position',pp1a(indr,:));
    SS.PlotUSDR(maxFR,[],sort_bound,0)
    set(gca,'XTick',[],'TickDir','out','TickLength',[0.005 0.025],'Box','off')
    xlabel ''
    
    %% Short time scale ASDR
%     axes('units','centimeters','position',pp1b(f,:)); 
%     set(gca,'XTick',[],'TickDir','out','TickLength',[0.005 0.025],'Box','off')
%     xlabel ''
    
    %% UNFR
    axes(aclamp);
    hold on
    plot([tTargOn/tdev tTargOff/tdev],[targetVal targetVal],'k:','LineWidth',1)
    stairs(SS.asdr.bin/tdev,SS.asdr.asdr/max(SS.unit),'color',col(ind,:),'LineWidth',2)
    xlim(tbound/tdev)
    ylim([0 7])
    ylabel('Hz/Unit')
    set(gca,'TickDir','out','TickLength',[0.005 0.025],'Xtick',[])
    
    % DUTY CYCLE
    axes(aduty);
    hold on
    t = SS.asdr.bin/tdev
    stairs(SS.asdr.bin/tdev,duty,'color',col(ind,:),'LineWidth',2)  
    ylim([0 maxduty+maxduty/2])
    xlim(tbound/tdev)
    ylabel('Duty Cycle')
    set(gca,'TickDir','out','TickLength',[0.005 0.025],'Box','off')
    xlabel('time (hrs)')
    
    drawnow
    
end
%%
axes(aduty);
plot([tTargOn/tdev tTargOff/tdev],[maxduty maxduty],'k:')

%%
% Export figure

% export_fig 'cl-12hr-bb-summary' -pdf -transparent
