% 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

% 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';
};
% fids = {
%     'multi-u_dead-time_2012-04-19_hoss-r3_tommy_358_01_nodrug';
%     'multi-u_dead-time_2012-04-19_hoss-r3_tommy_358_01_ap5';
%     };

close all
figure


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]);
    xlimit = [0 2500];
    
    %% Decode dig file
    b = dec2bin(d.state,32);
    u = bin2dec(b(:,1:29))/10e6;
    c = bin2dec(b(:,30:31));
    t = bin2dec(b(:,32));
    
    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
    
    %% Calculate ASDR for whole record
    ss.ASDR(0.1,[],[],[],[],0);
    
    %% For each unique u in uu find the average firing rate
    
    % num detected units
    nu = length(unique(ss.unit));
    fru = [];
    for i = 1:length(uu)
        fru = [fru mean(ss.asdr.asdr(ss.asdr.bin >= tu_on(i) & ss.asdr.bin <= tu_off(i)))/nu];
    end
    
    % fr_nom = mean(ss.asdr.asdr(ss.asdr.bin >= 0 & ss.asdr.bin <= min(tu_on)))/nu;
    fr_nom = mean(ss.asdr.asdr(ss.asdr.bin >= max(tu_on)+120))/nu;
    
    uFR = [];
    sFR = [];
    db = 0.1;
    bins = 0:db:0.9;
    
    % bin results
    for b = bins
        uFR = [uFR,mean(fru(uu>= b & uu < b + db))];
        sFR = [sFR,std(fru(uu>= b & uu < b + db))];
    end
    
    %% Plot result
    subplot(1,length(fids),f)
    semilogy(uu, fru, 'bo')
    hold on
    
    %% 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 + 20;
    
    % Calculate ASDR for whole record
    ss.ASDR(0.1,[],[],[],[],0);
    
    %% For each unique u in uu find the average firing rate
    
    % num detected units
    nu = length(unique(ss.unit));
    fru = [];
    for i = 1:length(uu)
        fru = [fru mean(ss.asdr.asdr(ss.asdr.bin >= tu_on(i) & ss.asdr.bin <= tu_off(i)))/nu];
    end
    [uuu ind] = sort(uu);
    fru = fru(ind);
    
    % fr_nom = mean(ss.asdr.asdr(ss.asdr.bin >= 0 & ss.asdr.bin <= tu_on(1)))/nu;
    fr_nom = mean(ss.asdr.asdr(ss.asdr.bin >= 0 & ss.asdr.bin < 600))/nu;
    
    uFR = [];
    sFR = [];
    db = 0.1;
    bins = 0:db:0.9;
    
    % bin results
    for b = bins
        uFR = [uFR,mean(fru(uuu>= b & uuu < b + db))];
        sFR = [sFR,std(fru(uuu>= b & uuu < b + db))];
    end
    
    %% Plot result
    subplot(1,length(fids),f)
    semilogy(uuu,fru, 'o','color',([255 205 42])/255)
    semilogy([0 1],[fr_nom fr_nom],'k-.')
    set(gca,'Fontsize',12,'TickDir','Out','TickLength',[0.015 0.025],'PlotBoxAspectRatio',[0.5 1 1]);
    axis([0 1 0.003 30])
    box off
    xlabel('u','Fontsize',12)
    ylabel('Hz/unit','Fontsize',12)
    
    
end

%% 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';
};

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]);
    xlimit = [0 2500];
    
    %% Decode dig file
    b = dec2bin(d.state,32);
    u = bin2dec(b(:,1:29))/10e6;
    c = bin2dec(b(:,30:31));
    t = bin2dec(b(:,32));
    
    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
    
    %% Calculate ASDR for whole record
    ss.ASDR(0.1,[],[],[],[],0);
    
    %% For each unique u in uu find the average firing rate
    
    % num detected units
    nu = length(unique(ss.unit));
    fru = [];
    for i = 1:length(uu)
        fru = [fru mean(ss.asdr.asdr(ss.asdr.bin >= tu_on(i) & ss.asdr.bin <= tu_off(i)))/nu];
    end
    
    % fr_nom = mean(ss.asdr.asdr(ss.asdr.bin >= 0 & ss.asdr.bin <= min(tu_on)))/nu;
    fr_nom = mean(ss.asdr.asdr(ss.asdr.bin >= max(tu_on)+120))/nu;
    
    uFR = [];
    sFR = [];
    db = 0.1;
    bins = 0:db:0.9;
    
    % bin results
    for b = bins
        uFR = [uFR,mean(fru(uu>= b & uu < b + db))];
        sFR = [sFR,std(fru(uu>= b & uu < b + db))];
    end
    
    %% Plot result
    subplot(1,length(fids),f)
    semilogy(uu, fru, 'b^', 'MarkerFaceColor','b')
    hold on
    
    %% 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 + 20;
    
    % Calculate ASDR for whole record
    ss.ASDR(0.1,[],[],[],[],0);
    
    %% For each unique u in uu find the average firing rate
    
    % num detected units
    nu = length(unique(ss.unit));
    fru = [];
    for i = 1:length(uu)
        fru = [fru mean(ss.asdr.asdr(ss.asdr.bin >= tu_on(i) & ss.asdr.bin <= tu_off(i)))/nu];
    end
    [uuu ind] = sort(uu);
    fru = fru(ind);
    
    % fr_nom = mean(ss.asdr.asdr(ss.asdr.bin >= 0 & ss.asdr.bin <= tu_on(1)))/nu;
    fr_nom = mean(ss.asdr.asdr(ss.asdr.bin >= 0 & ss.asdr.bin < 600))/nu;
    
    uFR = [];
    sFR = [];
    db = 0.1;
    bins = 0:db:0.9;
    
    % bin results
    for b = bins
        uFR = [uFR,mean(fru(uuu>= b & uuu < b + db))];
        sFR = [sFR,std(fru(uuu>= b & uuu < b + db))];
    end
    
    %% Plot result
    subplot(1,length(fids),f)
    semilogy(uuu,fru, '^','color',([255 205 42])/255, 'MarkerFaceColor',([255 205 42])/255)
    semilogy([0 1],[fr_nom fr_nom],'k-')
    set(gca,'Fontsize',12,'TickDir','Out','TickLength',[0.015 0.025],'PlotBoxAspectRatio',[0.5 1 1]);
    axis([0 1 0.003 30])
    box off
    xlabel('u','Fontsize',12)
    ylabel('Hz/unit','Fontsize',12)
    
    
end

%%
export_fig('-pdf','-transparent','combined_results_00')
