{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Script for fish logistic-growth harvesting project in Section 3.5.2.\n",
    "\n",
    "#The Population Data: The data from Table 3.10, starting with 1978, which we will call year 0.\n",
    "\n",
    "udata = [72148, 73793, 74082, 92912, 82323, 59073, 59920, 48789, 70638, 67462, 68702, 61191, 49599, \n",
    "         46266, 34877, 28827, 21980, 17463, 18057, 22681, 20196, 25776, 23796, 19240, 16495, 12167, \n",
    "         21104, 18871, 21241, 22962];"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#And the harvest rates from Table 3.10:\n",
    "hd = [0.18847, 0.149741, 0.21921, 0.17678, 0.28203, 0.34528, 0.20655, 0.33819, 0.14724, 0.19757, \n",
    "      0.23154, 0.20860, 0.33565, 0.29534, 0.33185, 0.35039, 0.28270, 0.19928, 0.18781, 0.19357, \n",
    "      0.18953, 0.17011, 0.15660, 0.28179, 0.25287, 0.25542, 0.08103, 0.087397, 0.081952, 0.10518];"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Number of data points, and years 0 to 29.\n",
    "n = len(udata);\n",
    "times = [i for i in range(n)];"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#A plot of the population\n",
    "pdata = list(zip(times,udata))\n",
    "plt1 = scatter_plot(pdata,facecolor='blue')\n",
    "plt1.axes_labels(['time (years)','Biomass'])\n",
    "show(plt1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#A plot of the harvest rates:\n",
    "pdata2 = list(zip(times,hd))\n",
    "plt2 = scatter_plot(pdata2);\n",
    "plt2.axes_labels(['time (years)','Harvest Rate'])\n",
    "show(plt2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Fitting the Data to the ODE: Form finite difference approximations (u(k+1) - u(k))/1 from the population\n",
    "#data, for k = 0 to k = 28). This approximates u'(t) when t = k:\n",
    "udiffdata = [udata[k+1]-udata[k] for k in range(n-1)];"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Next substitute u = udata[k] into the harvested logistic ODE u'(t) = r*u(t)*(1-u(t)/K) - h(t)*u(t)\n",
    "#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.\n",
    "#This expression depends on r and K.\n",
    "var('r K');\n",
    "udiffdata2 = [r*udata[k]*(1-udata[k]/K)-hd[k]*udata[k] for k in range(n-1)];"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Form a sum of squares that depends on r and K\n",
    "SS = function('SS')(r,K);\n",
    "SS(r,K) = add((udiffdata[k]-udiffdata2[k])^2 for k in range(n-2));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#To find the optimal r and K, start with a contour plot:\n",
    "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()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Something near r = 0.3, K = 200000 might be a good initial guess at a minimizer for SS.\n",
    "#If we use these values for r and K (but you should find the true minimizer) we could compare the\n",
    "#predicted cod population to the data by following the suggestion of Modeling Exercise 5.2.3.\n",
    "#Let U be an array to hold our numerical estimates of the cod population:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rr = 0.3; KK = 200000; #Temporary values for r and K\n",
    "U = [0 for k in range(n)];\n",
    "U[0] = udata[0];\n",
    "for k in range(n-1):\n",
    "    U[k+1] = U[k] + rr*U[k]*(1-U[k]/KK) - hd[k]*U[k];"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Plot the resulting predicted population along with the data\n",
    "pdatapred = list(zip(times,U))\n",
    "plt3 = scatter_plot(pdatapred,facecolor='red')\n",
    "plt3.axes_labels(['time (years)','Biomass'])\n",
    "show(plt1+plt3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "SageMath 9.2",
   "language": "sage",
   "name": "sagemath"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}