% Plot figures
% Updated by YB 29-08-18
% Previously, charge the eXX.mat data and fill this section
%load('Experiment34_180718_temporal_GFP_8h_ind_1h.mat')

%% Re-write conditions
nproteins = 1;      % number of fluorescent proteins
ReferencenM = 0;    % concentration of reference 1st, 2nd,etc
ControlStrain = 0;  % control of autofluorescence 1st, 2nd,etc
Gain_exp = 60;      % Gain of the experiment 
MEFL = 1;           % 121 molecules of fluorescein/cell for GFP
Line = 1.5;
Cmap = colormap(lines);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Concentrations and strains
k = 1;      %plot from this strain name
s = length(evol.medios.strains); %strains
induction = length(evol.medios.strains(k).data);      %inductions
replica = size(evol.medios.strains(k).data(k).OD, 2);  %replicas

if (ControlStrain~=0 )
replica_ControlStrain = size(evol.medios.strains(ControlStrain).data(k).OD, 2);  %replicas
end

time_total_0 = evol.time_OD - evol.time_OD(1);   % minutes
time_step = time_total_0(2) - time_total_0(1);   % minutes
x_pre = time_total_0(evol.npre).*ones(evol.npre,1);  % Time of pre_induction
Table_estimation = table;  %init table for save data


%% Figures
if nproteins < 2
    hfig = zeros(1,5);
else
    hfig = zeros(1,8);
end

for f = 1:size(hfig,2)
    hfig(f) = figure;
end

%% Initialization Axes and arrays
% For each strain
for j = 1:s
    % For each induction
    for n = 1:induction         
        %Control strain of autofluorescence
            %OD_control =  mean(evol.medios.strains(ControlStrain).data(n).OD,2) - evol.medios.OD_medio;
        if (ControlStrain~=0 )
            OD_control =  evol.medios.strains(ControlStrain).data(n).OD_medio - evol.medios.OD_medio;
            %F_control =  mean(evol.medios.strains(ControlStrain).data(n).F,2);
            F_control =  evol.medios.strains(ControlStrain).data(n).F_medio;
            if nproteins > 1
                F2_control =  evol.medios.strains(ControlStrain).data(n).F2_medio;  
            end
        else
            OD_control =  zeros(length(evol.medios.strains(k).data(n).OD),1);
            F_control =  zeros(length(evol.medios.strains(k).data(n).F),1);
        end
               
        %Corrected Absorbance and Fluorescence (without LB or M9)
        OD_0 =  evol.medios.strains(k).data(n).OD - evol.medios.OD_medio*ones(1,size(evol.medios.strains(k).data(n).OD,2));
        %F_0 = evol.medios.strains(k).data(n).F - OD_0.*((F_control./OD_control)*ones(1,size(OD_0,2)));
        
        %Checking F of control strain
        if k == ControlStrain
            F_0 = evol.medios.strains(k).data(n).F - evol.medios.F_medio*ones(1,size(OD_0,2));
        else 
            F_0 = evol.medios.strains(k).data(n).F - F_control*ones(1,size(OD_0,2));
        end
        
        if nproteins > 1
           if k == ControlStrain
                F2_0 = evol.medios.strains(k).data(n).F2 - evol.medios.F_medio*ones(1,size(OD_0,2));
           else 
                F2_0 = evol.medios.strains(k).data(n).F2 - F2_control*ones(1,size(OD_0,2));
           end            
        end
        
        %Reference concentration at 0nM
        if (n == ReferencenM && n~=0 )
            ODref = OD_0;
            Fref = F_0;
            if nproteins > 1
               F2ref = F2_0;
            end   
        else
            ODref = zeros(size(OD_0));
            Fref = zeros(size(F_0));
            if nproteins > 1
               F2ref = zeros(size(F2_0));
            end 
            
        end  
        
        %F/OD
        FOD = (F_0*MEFL)./(OD_0*Gain_exp);
        FODref = FOD- (Fref*MEFL)./(ODref*Gain_exp);
        if nproteins > 1
            F2OD = (F2_0*MEFL./OD_0*Gain_exp);
            F2ODref = F2OD-(F2ref*MEFL./ODref*Gain_exp);
        end
        
        %Save data for estimation only nprotein=1
        if n~=ReferencenM
            if j~=ControlStrain
            Id_name = evol.medios.strains(j).Id;
            Induction_name = evol.medios.strains(j).data(n).Inducc;
            ODname = ['OD_',Id_name,'_',Induction_name];
            FODname = ['FOD_',Id_name,'_',Induction_name];
            aux=table(FOD(1:end-1,:), OD_0(1:end-1,:),'Variablenames',{FODname, ODname});
            Table_estimation =[Table_estimation, aux];
            end
        end
                       
        %Plots       
        for f = 1:length(hfig)
            switch f
                case 1 %OD
                    figure(hfig(f));
                    hs = subplot(s,1,j);
                    hShEr = ShadedErrorBar(time_total_0, mean(OD_0,2), std(OD_0, 0, 2),{'LineWidth',Line,'Color', ColormapSeq6(n)}, 1);
                    hold on;

                    %Plotting induction instant
                    title(evol.medios.strains(k).Id);
                    xlabel(hs,'t (min)'); grid on;
                    ylabel(hs,'OD_{600}'); grid on;
                    if n == induction
                        hShErLine=findobj(hs,'Type','line','LineWidth',1.5);
                        legend(flipud(hShErLine), {evol.medios.strains(k).data.Inducc}); %flip the order
                    end
                       
                        if (n == induction && j==s)
                            hch = get(gcf,'Children');
                            set(hch(2:2:end),'XLim',[time_total_0(1) time_total_0(end)]); %the axes start in number 2
                            ymax = max([hch(2:2:end).YLim]);
                            ymin = min([hch(2:2:end).YLim]);
                            set(hch(2:2:end),'YLim',[0.01 ymax]);       
                            set(hch(2:2:end), 'YScale','linear');
                            set(hch(2:2:end), 'YTick',[0.01 0.1 0.5 ymax]);
                            OD_vector = ymin:(ymax-ymin)/(length(x_pre)-1):ymax;
                            for p=2:2:length(hch)
                                set(hch(p),'NextPlot','add');
                                plot(hch(p),x_pre, OD_vector,'--k');
                            end
                        end

                case 2 %F
                    figure(hfig(f));
                    hs = subplot(s,1,j);
                    hShEr = ShadedErrorBar(time_total_0, mean(F_0,2), std(F_0, 0, 2),{'LineWidth',Line,'Color', ColormapSeq6(n)}, 1);
                    hold on;
                    
                    %Plotting induction instant
                    title(evol.medios.strains(k).Id);
                    xlabel(hs,'t (min)'); grid on;
                    ylabel(hs,'F_1 (f.u.)'); grid on;
                    if n == induction
                        hShErLine=findobj(hs,'Type','line','LineWidth',1.5);
                        legend(flipud(hShErLine), {evol.medios.strains(k).data.Inducc} )
                    end
                    
                    if (n == induction && j==s)
                            hch = get(gcf,'Children'); 
                            set(hch(2:2:end),'XLim',[time_total_0(1) time_total_0(end)]);
                            ymax = max([hch(2:2:end).YLim]);
                            ymin = min([hch(2:2:end).YLim]);
                            set(hch(2:2:end),'YLim',[ymin ymax]);
                            OD_vector = ymin:(ymax-ymin)/(length(x_pre)-1):ymax;
                            for p=2:2:length(hch)
                                set(hch(p),'NextPlot','add');
                                plot(hch(p),x_pre, OD_vector,'--k');
                            end
                    end
                    
                case 3 %F1/A
                    figure(hfig(f));
                    hs = subplot(s,1,j);
                    hShEr = ShadedErrorBar(time_total_0, mean(FOD,2), std(FOD, 0, 2),{'LineWidth',Line,'Color', ColormapSeq6(n)}, 1);
                    hold on;
                    
                    %Plotting induction instant
                    title(evol.medios.strains(k).Id);
                    xlabel('t (min)'); grid on;
                    ylabel('F_1/A (r.f.u.)'); grid on;
                    if n == induction
                        hShErLine=findobj(hs,'Type','line','LineWidth',1.5);
                        legend(flipud(hShErLine), {evol.medios.strains(k).data.Inducc} )
                    end
                    
                    if (n == induction && j==s)
                            hch = get(gcf,'Children'); 
                            set(hch(2:2:end),'XLim',[time_total_0(1) time_total_0(end)]);
                            ymax = max([hch(2:2:end).YLim]);
                            ymin = min([hch(2:2:end).YLim]);
                            set(hch(2:2:end),'YLim',[ymin ymax]);
                            OD_vector = ymin:(ymax-ymin)/(length(x_pre)-1):ymax;
                            for p=2:2:length(hch)
                                set(hch(p),'NextPlot','add');
                                plot(hch(p),x_pre, OD_vector,'--k');
                            end
                    end
                     
                case 4 %F1/A wrt reference
                    figure(hfig(f));
                    hs = subplot(s,1,j);
                    hShEr = ShadedErrorBar(time_total_0, mean(FODref,2), std(FODref, 0, 2),{'LineWidth',Line,'Color', ColormapSeq6(n)}, 1);
                    hold on;
                    
                    %Plotting induction instant
                    title(evol.medios.strains(k).Id);
                    xlabel('t (min)'); grid on;
                    ylabel('F_1/A - F_r/A_r (r.f.u.)'); grid on;
                    if n == induction
                        hShErLine=findobj(hs,'Type','line','LineWidth',1.5);
                        legend(flipud(hShErLine), {evol.medios.strains(k).data.Inducc} )
                    end
                    
                    if (n == induction && j==s)
                            hch = get(gcf,'Children'); 
                            set(hch(2:2:end),'XLim',[time_total_0(1) time_total_0(end)]);
                            ymax = max([hch(2:2:end).YLim]);
                            ymin = min([hch(2:2:end).YLim]);
                            set(hch(2:2:end),'YLim',[ymin ymax]);
                            OD_vector = ymin:(ymax-ymin)/(length(x_pre)-1):ymax;
                            for p=2:2:length(hch)
                                set(hch(p),'NextPlot','add');
                                plot(hch(p),x_pre, OD_vector,'--k');
                            end
                    end
                     
                case 5 %OD wrt induction levels
                    figure(hfig(f));
                    hs = subplot(1,induction,n);
                    hShEr = ShadedErrorBar(time_total_0, mean(OD_0,2), std(OD_0, 0, 2),{'LineWidth',Line,'Color',Cmap(j,:)}, 1); 
                    hold on;
                    
                    %Plotting induction instant                   
                    xlabel(hs,'t (min)'); grid on;
                    ylabel(hs,'OD_{600}'); grid on;
                    if n == induction
                        hShErLine=findobj(hs,'Type','line','LineWidth',1.5);
                        legend(flipud(hShErLine), {evol.medios.strains(1:k).Id} )
                    end
                    
                    if j==s
                            %hlg = legend(set_hShEr,evol.medios.strains(k).Id,'Orientation','Vertical','Location','Best');
                            hch = get(gcf,'Children');
                            title(['OD at ' evol.medios.strains(1).data(n).Inducc 'nM']);
                            set(hch(2:2:end),'XLim',[time_total_0(1) time_total_0(end)]);
                            ymax = max([hch(2:2:end).YLim]);
                            ymin = min([hch(2:2:end).YLim]);
                            set(hch(2:2:end),'YLim',[ymin ymax]);                           
                            OD_vector = ymin:(ymax-ymin)/(length(x_pre)-1):ymax;
                            for p = 2:2:length(hch)
                                set(hch(p),'NextPlot','add');
                                plot(hch(p),x_pre, OD_vector,'--k');
                            end
                    end
               
                case 6
                    figure(hfig(f));
                    hs = subplot(1,s,j);
                    hShEr = ShadedErrorBar(time_total_0, mean(F2_0,2), std(F2_0, 0, 2),{'LineWidth',Line,'Color', ColormapSeq6(n)}, 1);
                    hold on;

                    %Plotting induction instant
                    title(evol.medios.strains(k).Id);
                    hlg = legend(set_hShEr6, evol.medios.strains(k).data(1:induction).Inducc,...
                            'Orientation','Vertical','Location','Best');
                    set(hlg,'FontSize',11);
                    xlabel('t (min)'); grid on;
                    ylabel('F_2 (f.u.)'); grid on;
                    
                    if (n == induction && j==s)
                            hch = get(gcf,'Children'); 
                            set(hch(2:2:end),'XLim',[evol.time_OD(1) evol.time_OD(end)]);
                            ymax = max([hch(2:2:end).YLim]);
                            ymin = min([hch(2:2:end).YLim]);
                            set(hch(2:2:end),'YLim',[ymin ymax]);
                            OD_vector = ymin:ymax;
                            %x_pre_vector = x_pre.*(ones(size(OD_vector)));
                            for p=2:2:length(hch)
                                set(hch(p),'NextPlot','add');
                                plot(hch(p),x_pre, OD_vector,'--k');
                            end
                    end
                                        
                case 7
                    figure(hfig(f));
                    hs = subplot(1,s,j);
                    hShEr = ShadedErrorBar(time_total_0, mean(F2OD,2), std(F2OD, 0, 2),{'LineWidth',Line,'Color', ColormapSeq6(n)}, 1);
                    hold on;
                    
                    %Plotting induction instant
                    title(evol.medios.strains(k).Id);
                    hlg = legend(set_hShEr7, evol.medios.strains(k).data(1:induction).Inducc,...
                        'Orientation','Vertical','Location','Best');
                    set(hlg,'FontSize',11);
                    xlabel('t (min)'); grid on;
                    ylabel('F_2/A (r.f.u.)'); grid on;
                    
                    if (n == induction && j==s)
                            hch = get(gcf,'Children'); 
                            set(hch(2:2:end),'XLim',[evol.time_OD(1) evol.time_OD(end)]);
                            ymax = max([hch(2:2:end).YLim]);
                            ymin = min([hch(2:2:end).YLim]);
                            set(hch(2:2:end),'YLim',[ymin ymax]);
                            OD_vector = ymin:ymax;
                            %x_pre_vector = x_pre.*(ones(size(OD_vector)));
                            for p=2:2:length(hch)
                                set(hch(p),'NextPlot','add');
                                plot(hch(p),x_pre, OD_vector,'--k');
                            end
                    end                 
            end
        end        
    end
    k = k+1;  %next strain
end


%%
aux=table(time_total_0(1:end-1),'VariableNames',{'Time_min'});
Table_estimation =[aux, Table_estimation];
estimation = table2struct(Table_estimation,'ToScalar',true);
save 'estimation' estimation;
%   %% Save figure
% Ancho = 8;
% Alto = 10;
% set(gcf,'units','centimeters');
% set(gcf,'papersize',[Ancho Alto]);
% set(gcf,'paperposition',[0,0,Ancho,Alto]);
% print -dpdf -painters -r300 'E34_4.pdf'; %djpg