clear;
close all;
[data, text]=xlsread('file_name',4);

data=data(~isnan(data(:,1)),:);% take only rows where the element in the firstz column is a reakl number

ndata=size(data,2)-1;

for i=1:ndata%11 data sets (excl time)
    G{1,i}=data(:,1)/60;
    G{2,i}=data(:,1+i);%-data(:,2);
end


dt=(G{1,i}(2)-G{1,i}(1));

smoothness=7;

tmin=10;

P2=zeros(2,ndata);
edge2=zeros(1,ndata);
ctr=[1:12];
% ser=[2,3];
% cys=[4];
% ile=[7,14];
% SerIle=[8,15];
% SerAsp=[6,13];
% AspIle=[9,16];
datasets=[ctr];

for i=datasets

lG{1,i}=G{1,i};% times
lG{2,i}=log(G{2,i}); % log of data

Gi{1,i}=lG{1,i};
Gi{2,i}=smooth(lG{1,i},lG{2,i},smoothness,'lowess');



dG{1,i}=lG{1,i}(2:end)-(dt/2);
dG{2,i}=(diff(Gi{2,i})/dt);%derivative of interpolated curve divided by O.D. = growth rate

%dG{1,i}=Gi{1,i};
%dG{2,i}=TVRegDiff(lG{2,i},10,10,[],'scale','large')/dt;



[P,~, ExtP]=mspeaks(dG{1,i},dG{2,i},'ShowPlot', 'true', 'Style','Extline','Denoising','true');%,'OverSegmentationFilter',20000);%true to view plots
if isempty(P)
    disp(['Skipping data set ',num2str(i),'. No peaks found']);
    continue;
end
II=P(:,1)>tmin;
P=P(II,:);
ExtP1=ExtP(II,:);
ExtP=ExtP1;

[~,i1]=sort(P(:,2),'descend');
P=P(i1,:);
ExtP=ExtP(i1,:);

Pmax(1,i)=P(1,1);%time of max growth
Pmax(2,i)=P(1,2);%max growth rate
edgemax(i)=ExtP(1,1);%left edge of peak
mOD(i)=spline(G{1,i},G{2,i},Pmax(1,i));%OD at max growth rate

eP=P(P(:,1)<P(1,1),:);%peaks before max peak
eExtP=ExtP(P(:,1)<P(1,1),:);

err=10^(-1);
[~,I1]=min(abs( dG{1,i}-edgemax(i)   ));% the index closest to the one of the time of the left edge
[~,I2]=min(abs( dG{1,i}-Pmax(1,i)   ));% the index closest to the one of the time of the maximum
hPmax=Pmax(2,i)-dG{2,i}(I1); %height of max peak from left edge to top

if ~isempty(eP) && eP(1,2)-dG{2,i}(  dG{1,i}==eExtP(1,2) )  >=  0.3*hPmax;
P2(1,i)=eP(1,1);
P2(2,i)=eP(1,2);%highest peak before max peak
edge2(i)=eExtP(1,1);%left edge of highest peak before max peak
end


x=abs(dG{2,i}-0.25*hPmax-dG{2,i}(I1)   );
[~,I]=min(x(I1:I2));%index of the closest point that is 0.25 of the max peak (search between the left edge and the maximum)
time25(i)=dG{1,i}(I+I1-1);

if edge2(i)>0  %there is a real second peak
    i
    avggrowth1(i)=mean(dG{2,i}( dG{1,i}>=edge2(i) & dG{1,i}<=Pmax(1,i)) );
    avggrowth2(i)=mean(dG{2,i}( dG{1,i}>=edge2(i) & dG{1,i}<=P2(1,i)) );
    avggrowth3(i)=mean(dG{2,i}( dG{1,i}>=edge2(i) & dG{1,i}<=time25(i) ));
    area(i)=dt*trapz(dG{2,i}( dG{1,i}>=edge2(i) & dG{1,i}<=time25(i) ));
else %there is no second peak
    avggrowth1(i)=mean(dG{2,i}( dG{1,i}>=edgemax(i) & dG{1,i}<=Pmax(1,i)) );
    avggrowth2(i)=0;
    avggrowth3(i)=mean(dG{2,i}( dG{1,i}>=edgemax(i) & dG{1,i}<=time25(i) ));
    area(i)=dt*trapz(dG{2,i}( dG{1,i}>=edgemax(i) & dG{1,i}<=time25(i) ));
end

dGOD{1,i}=spline(G{1,i},G{2,i},dG{1,i});%the O.D. at the timepoints of the growth rate dG; %the growth rate is 1/G dG/dt
dGOD{2,i}=dG{2,i};


%X = sprintf('%s: max growth rate %f  at time %d and OD %f. Average growth rate: %f',text{i+1},M,maxtime,mOD,r);
%disp(X)

% figure(10)
% plot(G{1,i},G{2,i},'LineWidth',2);

%pause;
end


out={'Condition','Lag time','Maximal Growth rate','Time of maximal growth rate','OD at max growth rate','time of left edge of minor peak','time of left edge of max peak','hight of minor peak','time of minor peak','time 25','Average growth rate from lag to max','average growth from lag to minor peak','Average growth from lag to time of 1/4 max peak height', 'area','Time from lag to max growth rate'};

for i=datasets
out{i+1,1}=text{i+1};
if edge2(i)>0
out{i+1,2}=edge2(i);
else
out{i+1,2}=edgemax(i);
end
out{i+1,3}=Pmax(2,i);
out{i+1,4}=Pmax(1,i);
out{i+1,5}=mOD(i);
out{i+1,6}=edge2(i);
out{i+1,7}=edgemax(i);
out{i+1,8}=P2(2,i);
out{i+1,9}=P2(1,i);
out{i+1,10}=time25(i);
out{i+1,11}=avggrowth1(i);
out{i+1,12}=avggrowth2(i);
out{i+1,13}=avggrowth3(i);
out{i+1,14}=area(i);
if edge2(i)>0
out{i+1,15}=Pmax(1,i)-edge2(i);
else
out{i+1,15}=Pmax(1,i)-edgemax(i);
end

end
xlswrite('growth.xls',out);



G1=G(:,datasets);
lG1=lG(:,datasets);
Gi1=Gi(:,datasets);
dG1=dG(:,datasets);
dGOD1=dGOD(:,datasets);

figure(1)
clf;

sp1=subplot(2,1,1);
p1=plot(G1{:});
set(p1,'Marker', 'none','LineStyle','-','LineWidth',2);
% legend([p1(1),p1(2),p1(3), p1(4)]);
% set(sp1,'colororderindex',2);
ylabel('O.D.');
xlabel('time (sec)');

sp2=subplot(2,1,2);
p2=plot(Gi1{:},'LineWidth',2);
hold on
plot(lG1{:},'Marker', '.','LineStyle','none');
hold off;
ylabel('Ln(O.D.)');
%set(sp2,'ColorOrderIndex',1);
xlabel('time (sec)');
% legend([p2(1),p2(2),p2(3), p2(4)]);


figure(2)
sp3=subplot(2,1,1);

p3=plot(dG1{:},'LineWidth',2);%,Gi1{1,i});
% set(sp3,'ColorOrderIndex',1);
% mean(dGs{2,4}(10:25)) % get a mean average of


ylabel('growth rate (sec^{-1})');
xlabel('time');

sp4=subplot(2,1,2);
p4=plot(dGOD1{:},'LineWidth',2);
%set(p4,'Marker', 'o');
% hold;
% set(sp4,'ColorOrderIndex',1);
% p5=plot(dGODs{:});
% set(p5,'Linewidth', 2);
% hold off;
ylabel('growth rate (sec^{-1})');
xlabel('O.D.');

% 

% 
% 
% figure(3)
% clf;
% 
% G2=G(:,7:end);
% lG2=lG(:,7:end);
% Gi2=Gi(:,7:end);
% dG2=dG(:,7:end);
% dGOD2=dGOD(:,7:end);
% 
% sp12=subplot(2,1,1);
% p12=plot(G2{:},'LineWidth',2);
% set(p12,'Marker', 'none','LineStyle','-');
% % legend([p1(1),p1(2),p1(3), p1(4)]);
% % set(sp1,'colororderindex',2);
% ylabel('O.D.');
% xlabel('time (sec)');
% 
% sp22=subplot(2,1,2);
% p22=plot(Gi2{:},'LineWidth',2);
% hold on
% plot(lG2{:},'Marker', '.','LineStyle','none');
% hold off;
% ylabel('Ln(O.D.)');
% % set(sp2,'ColorOrderIndex',1);
% xlabel('time (sec)');
% % legend([p2(1),p2(2),p2(3), p2(4)]);
% 
% 
% figure(4)
% sp32=subplot(2,1,1);
% 
% p32=plot(dG2{:},'LineWidth',2);%,Gi2{1,i});
% % set(sp3,'ColorOrderIndex',1);
% % mean(dGs{2,4}(10:25)) % get a mean average of
% 
% 
% ylabel('growth rate (sec^{-1})');
% xlabel('time');
% 
% sp42=subplot(2,1,2);
% p42=plot(dGOD2{:},'LineWidth',2);
% %set(p4,'Marker', 'o');
% % hold;
% % set(sp4,'ColorOrderIndex',1);
% % p5=plot(dGODs{:});
% % set(p5,'Linewidth', 2);
% % hold off;
% ylabel('growth rate (sec^{-1})');
% xlabel('O.D.');