#Script for fish logistic-growth harvesting project in Section 3.5.2.

#The Population Data: The data from Table 3.10, starting with 1978, which we will call year 0.

udata = [72148, 73793, 74082, 92912, 82323, 59073, 59920, 48789, 70638, 67462, 68702, 61191, 49599, 46266, 34877, 28827,
21980, 17463, 18057, 22681, 20196, 25776, 23796, 19240, 16495, 12167, 21104, 18871, 21241, 22962];

#And the harvest rates from Table 3.10:
hd = [0.18847, 0.149741, 0.21921, 0.17678, 0.28203, 0.34528, 0.20655, 0.33819, 0.14724, 0.19757, 0.23154, 0.20860,
0.33565, 0.29534, 0.33185, 0.35039, 0.28270, 0.19928, 0.18781, 0.19357, 0.18953, 0.17011, 0.15660, 0.28179, 0.25287,
0.25542, 0.08103, 0.087397, 0.081952, 0.10518];

#Number of data points, and years 0 to 29.
n = len(udata);
times = [i for i in range(n)];

#A plot of the population
pdata = list(zip(times,udata))
plt1 = scatter_plot(pdata,facecolor='blue')
plt1.axes_labels(['time (years)','Biomass'])
show(plt1)

#A plot of the harvest rates:
pdata2 = list(zip(times,hd))
plt2 = scatter_plot(pdata2);
plt2.axes_labels(['time (years)','Harvest Rate'])
show(plt2)

#Fitting the Data to the ODE: Form finite difference approximations (u(k+1) - u(k))/1 from the population
#data, for k = 0 to k = 28). This approximates u'(t) when t = k:
udiffdata = [udata[k+1]-udata[k] for k in range(n-1)];

#Next substitute u = udata[k] into the harvested logistic ODE u'(t) = r*u(t)*(1-u(t)/K) - h(t)*u(t)
#for k = 0 to k = n-1, to approximate the right side of the ODE at time t = k for k = 0 to k = n-1.
#This expression depends on r and K.
var('r K');
udiffdata2 = [r*udata[k]*(1-udata[k]/K)-hd[k]*udata[k] for k in range(n-1)];

#Form a sum of squares that depends on r and K
SS = function('SS')(r,K);
SS(r,K) = add((udiffdata[k]-udiffdata2[k])^2 for k in range(n-2));

#To find the optimal r and K, start with a contour plot:
contour_plot(log(SS(r,K)),(r,0.1,0.5),(K,40000,300000),aspect_ratio=0.000002,cmap='winter',fill=False,contours=50).show()

#Something near r = 0.3, K = 200000 might be a good initial guess at a minimizer for SS.
#If we use these values for r and K (but you should find the true minimizer) we could compare the
#predicted cod population to the data by following the suggestion of Modeling Exercise 5.2.3.
#Let U be an array to hold our numerical estimates of the cod population:

rr = 0.3; KK = 200000; #Temporary values for r and K
U = [0 for k in range(n)];
U[0] = udata[0];
for k in range(n-1):
    U[k+1] = U[k] + rr*U[k]*(1-U[k]/KK) - hd[k]*U[k];

#Plot the resulting predicted population along with the data
pdatapred = list(zip(times,U))
plt3 = scatter_plot(pdatapred,facecolor='red')
plt3.axes_labels(['time (years)','Biomass'])
show(plt1+plt3)