% 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

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 + 30;

        %% 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;
        FR_nom(j,f) = fr_nom;


    end
end


%% Optimal smoothing to create non-parametric models
u_halo = 0:0.01:1;
u_chr2 = 0:0.01:1;

mFR_ChR2 = zeros(2,4,length(u_chr2));
mFR_Halo = zeros(2,4,length(u_halo));
span = 0.90;

for k = 1:2
    for i = 1:length(fids)
        %         num = 75;
        %         spans = linspace(0.25,0.99,num);
        %         sse = zeros(size(spans));
        %
        %
        %         cp = cvpartition(length(U_ChR2{k,i}),'k',10);
        %         for j=1:length(spans),
        %             f = @(train,test) norm(test(:,2) - mylowess(train,test(:,1),spans(j)))^2;
        %             sse(j) = sum(crossval(f,[U_ChR2{k,i};FR_ChR2{k,i}]','partition',cp));
        %         end
        %
        %         [~,minj] = min(sse);
        %         span = spans(minj);
        
        mFR_ChR2(k,i,:) = mylowess([U_ChR2{k,i};FR_ChR2{k,i}]',u_chr2,span);
        
        
        %         cp = cvpartition(length(U_Halo{k,i}),'k',10);
        %         for j=1:length(spans),
        %             f = @(train,test) norm(test(:,2) - mylowess(train,test(:,1),spans(j)))^2;
        %             sse(j) = sum(crossval(f,[U_Halo{k,i};FR_Halo{k,i}]','partition',cp));
        %         end
        %
        %         [~,minj] = min(sse);
        %         span = spans(minj);
        
        mFR_Halo(k,i,:) = mylowess([U_Halo{k,i};FR_Halo{k,i}]',u_halo,span);
        
    end
end

%% Find the 90% values
u_95_chr2 = zeros(2,4);
u_95_halo = zeros(2,4);

fr_95_chr2 = zeros(2,4);
fr_95_halo = zeros(2,4);

p_cut = 1.0;

for k = 1:2
    for i = 1:length(fids)
        
        j_95 = find(max(mFR_Halo(k,i,:)) - mFR_Halo(k,i,:) >= p_cut*range(mFR_Halo(k,i,:)),1,'first');
        u_95_halo(k,i) = u_halo(j_95);
        fr_95_halo(k,i) = mFR_Halo(k,i,j_95)';
        
        j_95 = find(max(mFR_ChR2(k,i,:)) - mFR_ChR2(k,i,:) <= (1-p_cut)*range(mFR_ChR2(k,i,:)),1,'first');
        u_95_chr2(k,i) = u_chr2(j_95);
        fr_95_chr2(k,i) = mFR_ChR2(k,i,j_95)';
        
    end
end


%% Summarize Results
close all
figure('units','centimeter','position',[40 5 30 20]);
o1 = [5 8 3 3];
o2 = [5 4 3 3];
do = [4 0 0 0];
col = [ [0 100 255];
    [255 205 41]]/255;
ylimit1 = [0 20];
ylimit2 = [0 0.5];

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',3);
    plot(u_chr2, squeeze(mFR_ChR2(1,i,:)),'k')
    %     plot([0 1],[FR_nom(1,i) FR_nom(1,i)],'k:')
%     line([u_95_chr2(1,i), u_95_chr2(1,i),0],[0,fr_95_chr2(1,i),fr_95_chr2(1,i)],'Color',[1 0 0])
    
        plot(U_ChR2{2,i},FR_ChR2{2,i},'^','col',col(1,:),'MarkerSize',3);
    plot(u_chr2, squeeze(mFR_ChR2(2,i,:)),'k')
    %     plot([0 1],[FR_nom(2,i) FR_nom(2,i)],'k:')
%     line([u_95_chr2(2,i), u_95_chr2(2,i),0],[0,fr_95_chr2(2,i),fr_95_chr2(2,i)],'Color',[1 0 0])
    
    if i == 1
        set(gca,'Ylim',ylimit1,'Tickdir','out','TickLength',[0.025 1],'FontSize',8)
    else
        set(gca,'Ylim',ylimit1,'Tickdir','out','TickLength',[0.025 1],'FontSize',8)
    end
    
    axes('units','centimeter','position',o2+(i-1)*do);
    hold on
    plot(U_Halo{1,i},FR_Halo{1,i},'o','col',col(2,:),'MarkerSize',3);
    plot(u_halo, squeeze(mFR_Halo(1,i,:)),'k')
    %     plot([0 1],[FR_nom(1,i) FR_nom(1,i)],'k:')
%     line([u_95_halo(1,i), u_95_halo(1,i),0],[0,fr_95_halo(1,i),fr_95_halo(1,i)],'Color',[1 0 0])
    
        plot(U_Halo{2,i},FR_Halo{2,i},'^','col',col(2,:),'MarkerSize',3);
    plot(u_halo, squeeze(mFR_Halo(2,i,:)),'k')
    %     plot([0 1],[FR_nom(2,i) FR_nom(2,i)],'k:')
%     line([u_95_halo(2,i), u_95_halo(2,i),0],[0,fr_95_halo(2,i),fr_95_halo(2,i)],'Color',[1 0 0])
    
    if i == 1
        set(gca,'Ylim',ylimit2,'Tickdir','out','TickLength',[0.025 1],'FontSize',8)
    elseif i == 4
        set(gca,'Ylim',[0 5],'Tickdir','out','TickLength',[0.025 1],'FontSize',8)
    else
        set(gca,'Ylim',ylimit2,'Tickdir','out','TickLength',[0.025 1],'FontSize',8)
    end
    
end

% % figure('units','centimeter','position',[40 5 30 20]);
% o1 = [5 8 2.5 2.5];
% o2 = [5 5 2.5 2.5];
% ylimit1 = [0 20];
% ylimit2 = [0 0.25];
% 
% for i = 1:length(fids)
%     
%     axes('units','centimeter','position',o1+(i-1)*do);
%     hold on
%     plot(U_ChR2{2,i},FR_ChR2{2,i},'^','col',col(1,:),'MarkerSize',3);
%     plot(u_chr2, squeeze(mFR_ChR2(2,i,:)),'k')
%     %     plot([0 1],[FR_nom(2,i) FR_nom(2,i)],'k:')
%     line([u_95_chr2(2,i), u_95_chr2(2,i),0],[0,fr_95_chr2(2,i),fr_95_chr2(2,i)],'Color',[1 0 0])
%     
%     if i == 1
%         set(gca,'Ylim',ylimit1,'Tickdir','out','TickLength',[0.025 1])
%     else
%         set(gca,'Ylim',ylimit1,'Tickdir','out','TickLength',[0.025 1])
%     end
%     
%     axes('units','centimeter','position',o2+(i-1)*do);
%     hold on
%     plot(U_Halo{2,i},FR_Halo{2,i},'^','col',col(2,:),'MarkerSize',3);
%     plot(u_halo, squeeze(mFR_Halo(2,i,:)),'k')
%     %     plot([0 1],[FR_nom(2,i) FR_nom(2,i)],'k:')
%     line([u_95_halo(2,i), u_95_halo(2,i),0],[0,fr_95_halo(2,i),fr_95_halo(2,i)],'Color',[1 0 0])
%     
%     if i == 1
%         set(gca,'Ylim',ylimit2,'Tickdir','out','TickLength',[0.025 1])
%     elseif i == 4
%         set(gca,'Ylim',[0 5],'Tickdir','out','TickLength',[0.025 1])
%     else
%         set(gca,'Ylim',ylimit2,'Tickdir','out','TickLength',[0.025 1])
%     end
%     
% end

export_fig -pdf -transparent dynamic_range_30sec_fits_01






