% Jon Newman
% Multi-parameter activity assay for cultures transfected with
% optical exciters and silencers.
%
% This script combines data from both cultures on which this experiment was
% conducted
%   - venus-20708, 2012.05.03, 21 div, 17 dpt
%   - venus-20709, 2012.04.16, 18 div, 17 dpt

% Plot parameters
h = 22; % mm
w1 = 17; % mm
wg = 1.6; % mm
w2 = 3.4;
s = 3;
wt = w1+wg+w2+s;

close all
figure('Units','centimeters','Position',[1 1 30 10])

ap = {};
as = {};

for f = 1:4
    ap{f} = axes('Units','centimeters','Position',[5+(f-1)*wt/10 5 w1/10 h/10]);
    set(gca,'Fontsize',12,'TickDir','Out','TickLength',[0.015 0.025],'XTick',[0],'YScale','Log');
    if f == 1
        set(gca,'YTick',[0.01,0.1,1.0,10])
    else
        set(gca,'YTick',[])
    end
    axis([0 0.5 0.003 30])
    box off
    hold on
    
    as{f} = axes('Units','centimeters','Position',[5+(f-1)*wt/10+(w1+wg)/10 5 w2/10 h/10]);
    set(gca,'Fontsize',12,'TickDir','Out','TickLength',[0.015 0.025],'XTick',[1.0],'YTick',[],'YScale','Log');
    axis([0.9 1.0 0.003 30])
    box off
    hold on
end

%%


for j = 1:2
    switch j
        
        case {1}
            % Venus Data
            base = 'J:/newmonster/2012_05_03_Venus-MultiClamp/venus-20708/';
            
            fids = {
                'multi-u_dead-time_2012-05-03_hoss-r3_venus_20708_01';
                'multi-u_dead-time_2012-05-03_hoss-r3_venus_20708_01_cnqx';
                'multi-u_dead-time_2012-05-04_hoss-r3_venus_20708_01_ap5';
                'multi-u_dead-time_2012-05-08_hoss-r3_venus_20708_01_bic';
                };
            
            load('venus_20708_bad_unit.mat')
            
            sym1 = '^';
            sym2 = ':';
        case {2}
            %% Tommy data
            base = 'J:/newmonster/2012_04_16_Tommy-MultiClamp/tommy-20709/';
            
            fids = {
                'multi-u_dead-time_2012-04-16_hoss-r3_tommy_20709_01';
                'multi-u_dead-time_2012-04-17_hoss-r3_tommy_20709_01_cnqx';
                'multi-u_dead-time_2012-04-18_hoss-r3_tommy_20709_01_ap5';
                'multi-u_dead-time_2012-04-18_hoss-r3_tommy_20709_01_bic';
                };
            
            load('tommy_20709_bad_unit.mat')
            
            sym1 = 'o';
            sym2 = '-';
    end
    
    for f = 1:length(fids)
        
        
        fid = fids{f};
        s = loadspk([base fid '.spk']);
        d = loaddig([base fid '.dig']);
        ss = SqueakySpk(fid,s.meta.fs_Hz,1,s);
        ss.RemoveChannel([1 8 33 57 64 35 24 60]);
        ss.RemoveUnit([0 bu]);
        spk = ss.ReturnClean();
        
        
        %% Decode dig file
        b = dec2bin(d.state,32);
        u = bin2dec(b(:,1:29))/10e6;
        c = bin2dec(b(:,30:31));
        t = bin2dec(b(:,32));
        
        %% ChR2
        uu = unique(u(c == 0));
        uu = uu(uu~=0);
        
        % Find start and end time for all u
        tu_on = [];
        tu_off = [];
        for i = 1:length(uu)
            tu_on = [tu_on d.time(find(u == uu(i),1, 'first'))];
            tu_off = [tu_off d.time(find(u == uu(i),1, 'last'))];
        end
        
        % num detected units
        nu = max(ss.unit);
        fru = [];
        for i = 1:length(uu)
            dt =  (tu_off(i)-30) - tu_on(i);
            fru = [fru numel(spk.time(spk.unit~=0 & spk.time >= tu_on(i)& spk.time <= tu_off(i)))/dt/nu];
        end
        
        % Plot result
        semilogy(ap{f},uu(uu<=0.5), fru(uu<=0.5), sym1,'MarkerSize',3,'color',([000 100 255])/255)
        semilogy(as{f},uu(uu>=0.9), fru(uu>=0.9), sym1,'MarkerSize',3,'color',([000 100 255])/255)
        
        
        %% Halo
        uu = u(u~=0 & c == 1);
        
        % Find start and end time for all u
        tu_on = d.time(ismember(u,uu));
        tu_off = tu_on + 60;
        
        %% For each unique u in uu find the average firing rate
        
        % num detected units
        nu = max(ss.unit);
        fru = [];
        for i = 1:length(uu)
            dt =  tu_off(i) - tu_on(i);
            fru = [fru numel(spk.time(spk.unit~=0 & spk.time >= tu_on(i)& spk.time <= tu_off(i)))/dt/nu];
        end
        
        % Plot result
        semilogy(ap{f},uu(uu<=0.5), fru(uu<=0.5), sym1,'MarkerSize',3,'color',([255 205 42])/255)
        semilogy(as{f},uu(uu>=0.9), fru(uu>=0.9), sym1,'MarkerSize',3,'color',([255 205 42])/255)
        
        %% Fist 10 minutes of spontaneous firing to get nominal firing rate
        fr_nom = numel(spk.time(spk.unit~=0 & spk.time <= 600))/600/nu;
        
        % Plot result
        semilogy(ap{f},[0 0.5],[fr_nom fr_nom], sym2, 'Color',[0 0 0])
        semilogy(as{f},[0.9 1],[fr_nom fr_nom], sym2, 'Color',[0 0 0])
        
        drawnow
        
    end
end

%%
% saveas(gcf,'combined_results_alldrug_60sec.fig')
% export_fig('-pdf','-transparent','combined_results_alldrug_60sec')
