% Jon Newman
% 2013 04 1
% Examine the integrated optical delivery for different control types

% Load conversion table
load('./2012_03_23_trans-stim-r3_blue-yellow/iLED2pow.mat');

%% PI controller, no drugs
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';
    };

% close all
figure('units','centimeter','position',[5 5 20 20]);
a1 = axes('units','centimeter','position',[5 5 5 5]);
a2 = axes('units','centimeter','position',[12 5 5 5]);

for f = 1:length(fids)
    
    fid = fids{f};
    load(fid);
    
    % Convert control signals into average light powers
    stim_freq = avg_perf.uctl1*10+10;
    stim_pw   = avg_perf.uctl1*0.005;
    stim_duty = stim_freq.*stim_pw;
    stim_pow  = stim_duty.*(polyval(iLED2pow.p_blu,avg_perf.uctl1*1.5));
    semilogy(a1,avg_perf.steps(stim_pow > 10e-5),stim_pow(stim_pow > 10e-5),'ko-','MarkerSize',3)
    hold(a1,'on')
    
    
    % Convert control signals into average light powers
    stim_duty = 1;
    stim_pow  = stim_duty.*(polyval(iLED2pow.p_yel,avg_perf.uctl2*1.0));
    
    semilogy(a2,avg_perf.steps(stim_pow > 10e-5),stim_pow(stim_pow > 10e-5),'ko-','MarkerSize',3)
    hold(a2,'on')
    
end

%% On-off controller without drugs
fids = {
    'J:/newmonster/2012_02_09_Lollipop_Long-Term-optoBB/cl-long-bb_01_2012-02-09_chr2_lollipop_10803';      % 0.75 Hz
    'J:/newmonster/2012_02_10_Lollipop_Long-Term-optoBB/cl-long-bb_01_2012-02-10_chr2_lollipop_10803_1Hz';  % 1 Hz
    'J:/newmonster/2012_02_27_Lollipop_Long-Term-optoBB/2012_02_27_Lollipop_Long-Term-optoBB_00';           % 2 Hz
    'J:/newmonster/2012_02_11_Lollipop_Long-Term-optoBB/cl-long-bb_01_2012-02-11_chr2_lollipop_10803';      % 3 Hz
    'J:/newmonster/2012_02_16_Lollipop_Long-Term-optoBB/cl-long-bb_01_2012-02-16_chr2_lollipop_10803';      % 4 Hz
    'J:/newmonster/2012_02_14_Lollipop_Long-Term-optoBB/cl-long-bb_01_2012-02-14_chr2_lollipop_10803';      % 5 Hz
    'J:/newmonster/2012_02_12_Lollipop_Long-Term-optoBB/cl-long-bb_01_2012-02-12_chr2_lollipop_10803';      % 6 Hz
    };


targ = [];
pow  = [];
for f = 1:length(fids)
    
    fid = fids{f};
    
    % Load simoc file
    sim = loadsimoc([fid '.simoc']);
    
    % target value
    iTargOn = find(sim.target > 0,1,'first');
    iTargOff = find(sim.target > 0,1,'last');
    targetVal = sim.target(iTargOn);
    if targetVal == 1.5 % Account for code screw up early on
        targetVal = 1.0;
    end
    
    targ = [targ targetVal];
    
    % Stim power
    tduty = sim.time(sim.output(1,:) > 0);
    c = histc(tduty, 0:300:sim.time(end));
    stim_freq = 10*c/300/100; % convert to stim freq (update rate is 100 Hz, so need to account for that)
    stim_duty = mean(stim_freq(stim_freq~=0))*0.005;
    stim_pow  = stim_duty.*(polyval(iLED2pow.p_blu,2.0));
    pow  = [pow stim_pow];
    
end

semilogy(a1,targ,pow,'bs-','MarkerSize',3)

%% On-off controller with CNQX 1
sim = loadsimoc('J:\newmonster\2012_03_05_Lollipop-10803_Long-Term-optoBB_CNQX/2012_03_05_Lollipop-10803_Long-Term-optoBB_CNQX_00_20uM-CNQX.simoc');
targetVal = unique(sim.target(sim.target~=0));

% Stim power
tduty = sim.time(sim.output(1,:) > 0);
c = histc(tduty, 0:300:sim.time(end));
stim_freq = 10*c/300/100; % convert to stim freq (update rate is 100 Hz, so need to account for that)
stim_duty = mean(stim_freq(stim_freq~=0))*0.005;
stim_pow  = stim_duty.*(polyval(iLED2pow.p_blu,2.0));
semilogy(a1,targetVal,stim_pow,'r^','MarkerSize',3)

%% On-off controller with CNQX 2
bang = loadbang('Z:\jon\2013_02_26_Yam-22565_24-hr-cnqx-clamp\2013_02_26_Yam-22565_40-um-cnqx_01.bang');
targetVal = unique(bang.targ(bang.targ~=0));

tduty = bang.time([diff(bang.out1) 0] > 0);
c = histc(tduty, 0:300:bang.time(end));
stim_freq = c/300; % convert to duty cycle
stim_duty = mean(stim_freq(stim_freq~=0))*0.005;
stim_pow  = stim_duty.*(polyval(iLED2pow.p_blu,2.0));
semilogy(a1,targetVal,stim_pow,'g^','MarkerSize',3)


%% On-off controller with AP5 1
bang = loadbang('Z:\jon\2013_03_11_Yam-22564_24-hr-ap5-clamp\2013_03_11_Yam-22564_50um-ap5_01.bang');
targetVal = unique(bang.targ(bang.targ~=0));

tduty = bang.time([diff(bang.out1) 0] > 0);
c = histc(tduty, 0:300:bang.time(end));
stim_freq = c/300; % convert to duty cycle
stim_duty = mean(stim_freq(stim_freq~=0))*0.005;
stim_pow  = stim_duty.*(polyval(iLED2pow.p_blu,2.0));
semilogy(a1,targetVal,stim_pow,'yv','MarkerSize',3)

%% On-off controller with AP5 2
bang = loadbang('Z:\jon\2013_03_22_Yam-22565_24-hr-ap5-clamp\2013_03_22_Yam-22565_50um-ap5_01.bang');
targetVal = unique(bang.targ(bang.targ~=0));

tduty = bang.time([diff(bang.out1) 0] > 0);
c = histc(tduty, 0:300:bang.time(end));
stim_freq = c/300; % convert to duty cycle
stim_duty = mean(stim_freq(stim_freq~=0))*0.005;
stim_pow  = stim_duty.*(polyval(iLED2pow.p_blu,2.0));
semilogy(a1,targetVal,stim_pow,'mv','MarkerSize',3)

%% Touch up and save
set(a1,'TickDir','out','Ylim',[0.0005 2],'Xlim',[0 13])
set(a2,'TickDir','out','Ylim',[0.05 15],'Xlim',[0 13])

%%
saveas(gcf,'light_delivery_summary_02')
% export_fig -pdf -transparent light_delivery_summary_02_with-in-vivo_multi-scale.pdf



