% 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
% 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;
s   = 0.35; % verticle spacing
s2  = 0.15;

f1 = figure('units','centimeter','position',[40 5 30 40]);
a3 = axes('units','centimeter','position',[10 10 w h3]);
hold on;
a2 = axes('units','centimeter','position',[10 10+h3+s w h2]);
hold on;
tv = 10+h3+h2+2*s;
numsp = length(fids); % number of subplots
percent = 0;

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
    tTargOn = sim.time(iTargOn);
    tTargOff = sim.time(iTargOff);
    
    
    
    %% bang-bang input signal (duty cycle)
    
    maxduty = stimWidthSec*stimHz;
    tduty = sim.time([diff(sim.output(1,:)) 0] > 0)- t_start(1) + 3600*tdrug1 + 3600*tpre + 3600*tgap;
    c = histc(tduty, ss.asdr.bin);
    duty = c/dt; % convert to duty cycle
    
    % Plot Duty cycle
    b = duty~=0;
    plot(a3,SS.asdr.bin(b)/tdev,duty(b),'color',col(f,:))
    
    % Plot firing rates
    plot(a2, SS.asdr.bin/tdev,SS.asdr.asdr/max(SS.unit),'color',col(f,:))
    plot(a2, [tTargOn/tdev tTargOff/tdev],[targetVal targetVal],'k:')
    xlim(a2, tbound/tdev)
    
    
    %% USDR
    ax = axes('units','centimeter','position',[10 tv w h1x]);
    tv = tv + s2+h1x;
    SS.PlotUSDR(maxFR,[],sort_bound,0)
    set(gca,'XTick',[],'TickDir','out','TickLength',[0.005 0.025],'Box','off')
    xlabel ''
    ylabel ''
    
    percent = percent + 1/6;
    disp([percent ' % complete.'])
    drawnow
    
end

% Final plot additions
plot(a3, [0 tbound(2)/tdev],[maxduty maxduty],'k:')
set(a3,'Ylim',[0 maxduty+maxduty/4],'Xlim',tbound/tdev,'TickDir','out','TickLength',[0.005 0.025],'Box','off')
set(a2,'TickDir','out','TickLength',[0.005 0.025],'Xtick',[],'Ylim',[0 10])

%% Export figures

export_fig 'cl-long-bb_lollipop_10803_summary_02' -pdf -transparent
export_fig 'cl-long-bb_lollipop_10803_summary_02' -r400 -transparent