
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Fitting the logistic growth model to 4 sets of real data
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(1)

% linewidth for plots
lw=2;

% optional arguments for the optimization routines
%options=optimset('TolFun',1e-10,'TolX',1e-10,'MaxIter',1e10);
options=optimset;

% define the logistic growth 
% p(1) = r0 = growth rate
% p(2) = P0 = inital population
% p(3) = M  = carying capacity (a/b)
Pop = inline('p(3)./(1+exp(-p(1)*x)*(p(3)-p(2))/p(2))','p','x');

%%%%%%%%%%%%%%%%%%%%%%%% yeast data %%%%%%%%%%%%%%%%%%%%%%%%%%%%
% load data and store it in the X (time) and Y (population) vectors
data=load('yeast.dat');X=data(:,1);Y=data(:,2);

% define inital guess for the paramaters
r0=0.4; P0=10; M0=600; param0=[r0,P0,M0];

% apply a nonlinear square fit and store paramaters in pyeast
pyeast=lsqcurvefit(Pop,param0,X,Y,[],[],options)

% plot data and fitted solution
subplot(2,2,1)
t0=0; t1=max(X); T=t0:(t1-t0)/1000:t1;
plot(X,Y,'o',T,Pop(pyeast,T),'LineWidth',lw)
title('Yeast');
axis tight; set(gca,'LineWidth',1.5);

%%%%%%%%%%%%%%%%%%%%%%%% lupine data %%%%%%%%%%%%%%%%%%%%%%%%%%%%
data=load('lupine.dat');X=data(:,1);Y=data(:,2);
r0=0.4; P0=10; M0=200; param0=[r0,P0,M0];
plupine=lsqcurvefit(Pop,param0,X,Y,[],[],options)
subplot(2,2,2)
t0=0; t1=max(X); T=t0:(t1-t0)/1000:t1;
plot(X,Y,'o',T,Pop(plupine,T),'LineWidth',lw)
title('Lupine')
axis tight; set(gca,'LineWidth',1.5);

%%%%%%%%%%%%%%%%%%%%%%%% bacteria data %%%%%%%%%%%%%%%%%%%%%%%%%%%%
data=load('bacteria.dat');X=data(:,1);Y=data(:,2);
r0=0.4; P0=10; M0=60; param0=[r0,P0,M0];
pbacteria=lsqcurvefit(Pop,param0,X,Y,[],[],options)
subplot(2,2,3)
t0=0; t1=max(X); T=t0:(t1-t0)/1000:t1;
plot(X,Y,'o',T,Pop(pbacteria,T),'LineWidth',lw)
title('Bacteria')
axis tight; set(gca,'LineWidth',1.5);

%%%%%%%%%%%%%%%%%%%%%%%% bombay data %%%%%%%%%%%%%%%%%%%%%%%%%%%%
data=load('bombay.dat');X=data(:,1);Y=data(:,2);
r0=0.4; P0=10; M0=10000; param0=[r0,P0,M0];
pbombay=lsqcurvefit(Pop,param0,X,Y,[],[],options)
subplot(2,2,4)
t0=0; t1=max(X); T=t0:(t1-t0)/1000:t1;
plot(X,Y,'o',T,Pop(pbombay,T),'LineWidth',lw)
title('Bombay')
axis tight; set(gca,'LineWidth',1.5);

% save plot in a JPG file:
%print -djpeg  bombay_yeast_lupine_bacteria.jpg


