% 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

%% Data storage
FR_nom = zeros(2,4); % spont firing rates
FR_nom_units = cell(2,4); % spont firing rates

U_ChR2 = cell(2,4);  % Control signals
U_Halo = cell(2,4);

FR_ChR2 = cell(2,4); % Evoked firing rates
FR_Halo = cell(2,4);

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)-30)))/dt/nu];
        end


        U_ChR2{j,f} = uu';
        FR_ChR2{j,f} = fru;


        %% 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

        U_Halo{j,f} = uu';
        FR_Halo{j,f} = fru;

        %% Fist 10 minutes of spontaneous firing to get nominal firing rate
        fr_nom = numel(spk.time(spk.unit~=0 & spk.time <= 600))/600/nu;
        units = unique(spk.unit(spk.unit~=0));
        fr = zeros(size(units));
        for u = 1:length(units)
            
            fr(u) = numel(spk.time(spk.unit== units(u) & spk.time <= 600))/600;
        end

        FR_nom(j,f) = fr_nom;
        FR_nom_units{j,f} = fr;


    end
end

%% Load LED power data
load('../2012_03_23_trans-stim-r3_blue-yellow/iLED2pow.mat');

%% Combine data into matracies (averaging over cultures)
U_H     = cell2mat(U_Halo');
FR_H    = cell2mat(FR_Halo');
U_C     = cell2mat(U_ChR2');
FR_C    = cell2mat(FR_ChR2');

%% Optimal smoothing to create non-parametric models
u_eval = 0:0.01:1;
frac = 0.95;

fit_FR_ChR2 = zeros(4,length(u_chr2));
p_fit_ChR2 = cell(4,1);
fit_95_ChR2 = zeros(4,4);
fit_FR_Halo = zeros(4,length(u_halo));
p_fit_Halo = cell(4,1);
fit_95_Halo = zeros(4,2);
d_range = zeros(4,2);

close all
for i = 1:size(U_H,1)
    
    start = [0.610958658746201 0.778802241824093 0.266471490779072];
    p_fit_ChR2{i} = create_chr2_fit(U_C(i,:),FR_C(i,:),1,start);
    fit_FR_ChR2(i,:) = p_fit_ChR2{i}(u_eval);
    u_sat = (-p_fit_ChR2{i}.C^-1)*log((1/frac -1)/p_fit_ChR2{i}.B) - 1;
    
    % U_sat, freq., width, 470 nm light power
    fit_95_ChR2(i,:) = [u_sat, 10 + 10*u_sat,5*u_sat ,polyval(iLED2pow.p_blu,u_sat)];
    
    d_range(i,1) = p_fit_ChR2{i}(1);
    
    start = -[0.610958658746201 0.778802241824093 0.266471490779072  0.0679927684700106];
    p_fit_Halo{i}= create_halo_fit(U_H(i,:),FR_H(i,:),1,start);
    fit_FR_Halo(i,:) = p_fit_Halo{i}(u_eval);
    fr_0 = p_fit_Halo{i}(0);
    u_sat = -p_fit_Halo{i}.C*log(((1 - frac)*(fr_0 - p_fit_Halo{i}.D))/p_fit_Halo{i}.A) + p_fit_Halo{i}.B;
    
    % U_sat, 590 nm light power
    fit_95_Halo(i,:) = [u_sat, polyval(iLED2pow.p_yel,u_sat)];
    
    d_range(i,2) = p_fit_Halo{i}(1);
    
end

%% Summarize Results
close all
figure('units','centimeter','position',[40 5 30 20]);
ms = 2;
o1 = 1.15*[5 8 1.2 1.5];
do = [6 0 0 0];

col = [ [0 100 255];
    [255 205 41]]/255;
ylimit = [0.005 30];
xlimit = [0 0.6];

for i = 1:length(fids)
    
    axes('units','centimeter','position',o1+(i-1)*do);
    hold on
    plot(U_ChR2{1,i},FR_ChR2{1,i},'o','col',col(1,:),'MarkerSize',ms);
    plot(U_ChR2{2,i},FR_ChR2{2,i},'^','col',col(1,:),'MarkerSize',ms);
    plot(u_eval, fit_FR_ChR2(i,:),'b-')
    
    plot(U_Halo{1,i},FR_Halo{1,i},'o','col',col(2,:),'MarkerSize',ms);
    plot(U_Halo{2,i},FR_Halo{2,i},'^','col',col(2,:),'MarkerSize',ms);
    plot(u_eval, fit_FR_Halo(i,:),'b-')
    
    plot([0 1],[FR_nom(1,i) FR_nom(1,i)],'r-')
    plot([0 1],[FR_nom(2,i) FR_nom(2,i)],'r:')
    
    if i == 1
        set(gca,'Ylim',ylimit,'Xlim',xlimit,'Tickdir','out','TickLength',[0.025 1],'FontSize',7,'YScale','log','YTick',[0.01 0.1 1 10],'XTick',[0 0.2 0.4 0.6])
    else
        set(gca,'Ylim',ylimit,'Xlim',xlimit,'Tickdir','out','TickLength',[0.025 1],'FontSize',7,'YScale','log','YTick',[0.01 0.1 1 10],'YTickLabel',{},'XTick',[0 0.2 0.4 0.6])
    end
    
end

o1 = 1.15*[5+1.3 8 0.2 1.5];
do = [6 0 0 0];
xlimit = [0.9 1.0];

for i = 1:length(fids)
    
    axes('units','centimeter','position',o1+(i-1)*do);
    hold on
    plot(U_ChR2{1,i},FR_ChR2{1,i},'o','col',col(1,:),'MarkerSize',ms);
    plot(U_ChR2{2,i},FR_ChR2{2,i},'^','col',col(1,:),'MarkerSize',ms);
    plot(u_eval, fit_FR_ChR2(i,:),'b-')
    
    plot(U_Halo{1,i},FR_Halo{1,i},'o','col',col(2,:),'MarkerSize',ms);
    plot(U_Halo{2,i},FR_Halo{2,i},'^','col',col(2,:),'MarkerSize',ms);
    plot(u_eval, fit_FR_Halo(i,:),'b-')
    
    plot([0 1],[FR_nom(1,i) FR_nom(1,i)],'r-')
    plot([0 1],[FR_nom(2,i) FR_nom(2,i)],'r:')
    
    set(gca,'Ylim',ylimit,'Xlim',xlimit,'Tickdir','out','TickLength',[0.025 1],'FontSize',7,'YScale','log','YTick',[],'XTick',[1.0])
    
end


saveas(gcf,'dynamic_range_60sec_fits_02')
export_fig -pdf -transparent dynamic_range_60sec_fits_02






