% Jon Newman

%% file id
fids = { 
    'J:/newmonster/2012_07_02_Xeko-1-minute-clamps/2012_07_02_Xeko-20709_one-min-clamp_07.mat';
    'J:/newmonster/2012_10_11-Herb-I5387-1-minute-clamps/2012_10_11_Herb-I5387_one-min-clamp_02.mat';
    'J:/newmonster/2012_10_22_Herb-I5386-1-minute-clamps/2012_10_22_Herb-I5386_one-min-clamp_10.mat';
    'Z:/jon/2013_02_14_Wingnut-22567-1-minute-clamps/instability-test/2013-02-14_wingnut_22567_one-minute-clamp_07.mat'; % Unstable network
    'Z:/jon/2013_02_26_Yam-22565_1-minute-clamp/2013-02-26_yam_22565_one-minute-clamp_06.mat';
    'Z:/jon/2012_11_21_Mandelbrot-21332-1-minute-clamps/2012_11_21_Mandelbrot_21332_one-min-clamp_12.mat';
    'Z:/jon/2013_04_14_Dojo-i5407_1-minute-clamps-with-blockers/2013_04_14_Dojo-i5407_one-min-clamp_nodrug_05.mat';
    };



sym1 = {'o','o','o','o','o','o'};


RT = nan(length(fids),13); % rise-time
TAU = nan(length(fids),13); % time-constant
ST = nan(length(fids),13); % settling time
JIT = nan(length(fids),13); % settling time
% ctl_2 = nan(length(fids),11);

dfr_r = 0.1;
dfr_s = 0.25;
sw = 2.5; %sec
st_time = 2; % second
loop_dt = 0.01;

for f = 1:length(fids)
    
    fid = fids{f};
    load(fid);
    
    % Determine rise-time measures
    for i = 1:length(avg_perf.steps)
        
        targ = avg_perf.steps(i);
        
        b = ts.t{i} > 0;
        t = ts.t{i}(b);
        loop_dt = mean(diff(t));
        
        tru_fr = ts.fr{i}(b);
        sm_fr = smooth(ts.fr{i}(b),round(sw/loop_dt),'lowess');
        
        if size(sm_fr) ~= size(t)
            sm_fr = sm_fr';
        end
        
        
        if targ >= tru_fr(1) 
            rt_i = find(tru_fr - targ > -dfr_s & t > 0.05, 1 ,'first');
        else
            rt_i = find(tru_fr - targ < dfr_s & t > 0.05, 1 ,'first');
        end
        
        
      
        
        
        
        if ~isempty(rt_i)
            RT(f,i) = loop_dt*rt_i;
        end

        
        ind_in_bound = find(abs(sm_fr - targ) < dfr_s); % & t > 2);
        ind_out_bound = find(abs(sm_fr - targ) > dfr_s);
        
        if  isempty(ind_in_bound)
            continue;
        elseif size(ind_in_bound,2) > 1
            ind_in_bound = ind_in_bound';
            ind_out_bound = ind_out_bound';
        end
        
        ind_enter_bound = [ind_in_bound(1);  ind_in_bound(find(diff(ind_in_bound) > 1)+1)];
        ind_exit_bound = [ind_out_bound(1); ind_out_bound(find(diff(ind_out_bound) > 1)+1)];
        
        
        if ind_exit_bound(1) < ind_enter_bound(1)
            ind_exit_bound(1) = [];
        end
        if length(ind_exit_bound) ~= length(ind_enter_bound)
            ind_exit_bound = [ind_exit_bound; length(sm_fr)];
        end
        
        t_in_bound = ind_exit_bound-ind_enter_bound;
        ind_settled = ind_enter_bound(find(t_in_bound > st_time/loop_dt,1,'first'));
        
        if ~isempty(ind_settled)
            ST(f,i) = loop_dt*ind_settled;
            JIT(f,i) = avg_perf.sfr(i);

        end
        
        close all
        figure
        hold on
        plot(tru_fr,'y')
        plot(sm_fr,'k')
        plot(ind_settled,sm_fr(ind_settled),'r*')
        plot(rt_i,sm_fr(rt_i),'g*')
      
        
    end
    
end

ST(JIT > 0.5) = nan;

JIT(JIT > 0.5) = nan; % remove the case where the there was a control failure


%%

close all

col = [ [0 100 255];
    [255 205 41]]/255;

col2 = [
    [220 000 000];
    [255 102 000];
    [255 205 041];
    [000 170 000];
    [000 170 212];
    [000 100 255];
    [085 000 212];
    [102 102 102]
    ]/255;

figure('units','centimeter','position',[40 5 30 20]);

ax1 = axes('units','centimeter','position',[5 10 6 3]);
hold on;
set(gca,'TickDir','out','FontName','Arial','FontSize',10,'TickLength',[0.025 0.025],'Xlim',[0 10],'Ylim',[0 25]);
ylabel('Settling time (sec)')
xlabel('Target rate (Hz/unit)')


for i = 1:size(ST,1)
    x = 0:10;
    plot(x,ST(i,1:11),'o','MarkerSize',4,'Color',col2(i,:),'MarkerSize',3,'MarkerFaceColor',[1 1 1]);
    ind_out = find(ST(i,1:11) > 25);
    if ~isempty(ind_out)
       plot(x(ind_out),24*ones(size(ind_out)), 'r^')
        
    end
    
end

plot(0:12,nanmean(ST),'k-','Linewidth',2)
plot(0:12,nanmedian(ST),'k:','Linewidth',2)

% Median and CI's
conf = .95;
a = 1-conf;
b1 = norminv(1-a/2);
b2 = b1/sqrt(2);
rfac = (1/(norminv(.75)-norminv(.25))) * sqrt(pi/2) * (b1+b2)/2;

dat = reshape(ST(~isnan(ST)),1,[]);
n = size(dat,2);
q = quantile(dat,[0.25, 0.5, 0.75]);
q(2)
ci_u = q(2) + rfac*(q(3)-q(1))/sqrt(n)
ci_l = q(2) - rfac*(q(3)-q(1))/sqrt(n)



% ax2 = axes('units','centimeter','position',[10 10 3 3]);
% hold on;
% set(gca,'TickDir','out','FontName','Arial','FontSize',10,'TickLength',[0.025 0.025],'Xlim',[0 10],'Ylim',[0 0.2]);
% ylabel('SD FR (Hz/unit)')
% xlabel('Target rate (Hz/unit)')
% 
% for i = 1:size(ST,1)
%     
%     plot(0:10,JIT(i,1:11),'o','MarkerSize',4,'Color',col2(i,:),'MarkerSize',3,'MarkerFaceColor',[1 1 1])
% end
% 
% plot(0:12,nanmean(JIT),'k-','Linewidth',2)


% ax3 = axes('units','centimeter','position',[15 10 3 3]);
% ax2 = axes('units','centimeter','position',[10 10 3 3]);
% % hold on;
% 
% 
% for i = 1:size(RT,1)
%     
%     plot(0:10,RT(i,1:11),'o','MarkerSize',4,'Color',col2(i,:),'MarkerSize',3,'MarkerFaceColor',[1 1 1])
%     hold on
% end
% 
% plot(0:12,nanmean(RT),'k-','Linewidth',2)
% set(gca,'TickDir','out','FontName','Arial','FontSize',10,'TickLength',[0.025 0.025],'Xlim',[0 10],'Ylim',[0.01 10]);
% ylabel('SD FR (Hz/unit)')
xlabel('Target rate (Hz/unit)')


%%
saveas(gca,'1-minute-clamp_settling-time_and_jitter_02.fig')
export_fig -pdf -transparent 1-minute-clamp_settling-time_and_jitter_02.pdf

