%% Fish Harvesting Revisited
%% Kurt Bryan and SIMIODE
% Script for fish logistic-growth harvesting project in Section 3.5.2.
%
% The Data: The data from Table 3.10, starting with 1978, which we will call year 0. Note Matlab indexes arrays starting with index 1, however.

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];

% 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];

% Plot of the population data, indexed as years 0 to 29:

plot([0:29],udata)

% And a plot of the harvest rates

plot([0:29],hd)

% Fitting the Data to the ODE: Form finite difference approximations (u(k+1) - u(k))/1 from the population data, for k = 1 to k = 29). This approximates u'(t) when t = k:

udiffdata = udata(2:30)-udata(1:29);

% 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.

syms r;
syms K;
u1 = udata(1:29);
udiffdata2 = r*u1.*(1-u1/K)-hd(1:29).*u1;

% Form a sum of squares that depends on r and K

SS(r,K) = sum((udiffdata-udiffdata2).^2);

% Finding the optimal r and K: Start with a plot of SS(r, K), or log(SS).

fsurf(log(SS), [0.1 0.5 40000 300000])

% Or perhaps a contour plot, with contours in steps of 0.05 from 21.5 to 24:

fcontour(log(SS), [0.1 0.5 40000 300000],'LevelList',[21.5:0.05:24])

% 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:

U = zeros(1,30);
U(1) = udata(1);

% Now use the harvested logistic ODE with our estimates for r and K to march the predicted cod biomass out in time (this is essentially Euler's method):

r = 0.3;
K = 200000;
for k=1:29
    U(k+1) = U(k) + r*U(k)*(1-U(k)/K) - hd(k)*U(k);
end

% Plot the resulting predicted population along with the data

plot([0:29],udata)
hold on
plot([0:29],U,'-r')
hold off

%
%