#Competing Yeast Species

#A worksheet to facilitate analyzing the competing yeast species project in Section
#7.5.3 of the text. The goal is to use the data of Table 7.2 to estimate the
#parameters r1, K1, r2, K2, a, and b in equations (7.3)-(7.4).

#Load in the data. Missing population data is indicated with a -1.
load('comp_yeast_data.sage');

#Let's form amalgamated data sets for runs 1 and 2, for the Saccharomyces population
#alone (column 2), the Schizosaccharomyces population alone (column 4), the
#Saccharomyces population in competition (column 3), and the Schizosaccharomyces
#population in competition (column 5).

N0 = N1 + N2; #Total number of times/data points.
times = times1+times2; #Concatenate all times
sacc = sacc1+sacc2; #Concatenate Sacch data
saccind = list(zip(times,sacc)); #Arrange in (time, population) pairs.
schiz = schiz1+schiz2; #Concatenate Schiz data
schizind = list(zip(times,schiz)); #Arrange in (time, population) pairs.
saccm = saccmix1+saccmix2; #Concatenate Sacch data
saccmix = list(zip(times,saccm)); #Arrange in (time, population) pairs.
schizm = schizmix1+schizmix2; #Concatenate Schiz data
schizmix = list(zip(times,schizm)); #Arrange in (time, population) pairs.

#We can plot the data. For example, for the saccharomyces population growing in
#isolation. Also, limit view to 0 to 60 in time, 0 to 15 in population.
list_plot(saccind).show(xmin=0,xmax=60,ymin=0,ymax=15);

#A similar plot can be constructed for the schizosaccharomyces population in
#isolation, or either species in competition.


#Estimate Growth Parameters for Saccharomyces: We can estimate the parameters
#r1 and K1 for the logistic growth model for the saccharomyces population
#using the corresponding data.

#The solution to the logistic ODE u' = r1*u*(1-u/K1) with initial data u(t0) = u0 is
var('t r1 K1 u0 t0'); #Declare t as symbolic variable, t0,u0,r1,K1 too.
usol(t) = K1*u0/(exp(-r1*(t-t0))*(K1-u0)+u0);

#The initial data could be taken as any (time, population) data point, but
#we'll use the first:
t0 = times[0];
u0 = sacc[0];

#The solution with this initial data is
u1sol(t) =  K1*u0/(exp(-r1*(t-t0))*(K1-u0)+u0);

#A plot or even a glance at the data suggests that K1 is around 12 to 14. To
#estimate K1 and r1 more carefully,  we can use the "find_fit" command, that performs
#a least squares fit of u1sol(t) to the data by adjusting r1 and K1.

#But first, we must filter out the "no data" points by excluding those with population = -1
saccpos =  [[times[i],sacc[i]] for i in range(N0) if sacc[i]>0]

#Now fit the model
sol = find_fit(saccpos,u1sol,parameters=[r1,K1],initial_guess=[0.5,14.0]) #Fit the model by adjusting r1 and K1
r1best = sol[0].rhs();
K1best = sol[1].rhs();
bestu1(t) = u1sol(r1=r1best,K1=K1best) #Define the best fit function of the form in "model"

#Plot the best fit u1(t) and the data.
p1 = list_plot(saccind);
p2 = plot(bestu1(t),(t,0,60),color='red');
(p1+p2).show(xmin=0,xmax=60,ymin=0,ymax=15);

#A similar procedure can be used to estimate r2 and K2 for the schizosaccharomyces data,
#using (for example) the data point t0 = 16, u0 = 1 as an "initial" data point.

#The Competition Parameters a and b: The ODEs that govern the competing species system are
var('a b r2 K2');
u1 = function('u1')(t);
u2 = function('u2')(t);
de1 = diff(u1,t) == r1*u1*(1-u1/K1 - a*u2/K1)
de2 = diff(u2,t) == r2*u2*(1-u2/K2 - b*u1/K2)

#A closed-form solution cannot be obtained, but a guess-and-plot approach based on
#solving these ODEs numerically can work, using the values for r1, K1, r2, K2 found above.