{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "63f5174f",
   "metadata": {},
   "source": [
    "# Modul 3 Kodeeksempler forskerforedrag"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "314c2bb5",
   "metadata": {},
   "source": [
    "## Opsætning"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f48abc98",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.integrate import solve_ivp\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib\n",
    "if not hasattr(matplotlib.RcParams, \"_get\"):\n",
    "    matplotlib.RcParams._get = dict.get"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fc4f46c7",
   "metadata": {},
   "source": [
    "## Solve system 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2245ed90",
   "metadata": {},
   "outputs": [],
   "source": [
    "# --------------------------------------------------------------------\n",
    "# Tragedy of the commons generic ODE system - math version\n",
    "# --------------------------------------------------------------------\n",
    "# state vector y = (B,A)\n",
    "# B: biomass\n",
    "# A: actors\n",
    "# --- biological parameters ---\n",
    "# r: growth rate              \n",
    "# K: carrying capacity      \n",
    "# --- economical parameters ---\n",
    "# s: sector adaptation rate \n",
    "# p: catch price             \n",
    "# c: operating cost\n",
    "#\n",
    "# In this version parameter values are close to unity and not so\n",
    "# closely reflecting real world systems - it has 5 parameters.\n",
    "# By non-dimensionalization the system can further be reduced to\n",
    "# 3 essential parameters. Parameter/variable units are not spelled out.\n",
    "#\n",
    "# It has 2 parameter settings - one with an interior equilibrium, and one with\n",
    "# a boundary equilibrium (sector crash)\n",
    "# --------------------------------------------------------------------\n",
    "def ToC(t,y, r,K,p,c,s):\n",
    "    B,A = y\n",
    "    dBdt = r*B*(1-B/K) - A*B\n",
    "    dAdt = s*(p*B/c - 1)\n",
    "    if A<0:\n",
    "        dAdt = max(dAdt, 0)     # dAdt > 0 is OK if A<0, allow fleet to rebuild\n",
    "    return np.array([dBdt,dAdt])\n",
    "\n",
    "# ---- this parameter combination corresponds to an interior equilibrium ----\n",
    "r  = 1    # growth rate\n",
    "K  = 1    # carrying capacity\n",
    "p  = 1    # catch price  \n",
    "c  = 0.6  # operating cost\n",
    "s  = 1.5  # sector adaptation rate \n",
    "y0 = (K, 0) # initial condition\n",
    "\n",
    "# ---- this parameter combination is in sector crash regime; uncomment line below\n",
    "# c   = 1.4; y0  = (1.3*K, 1) \n",
    "\n",
    "tspan  = (0, 20)\n",
    "t_eval = np.linspace(tspan[0], tspan[1], 1000)\n",
    "args   = (r,K,p,c,s)\n",
    "sol1   = solve_ivp(ToC, tspan, y0, t_eval=t_eval, args=args)\n",
    "\n",
    "# -------- plotting --------\n",
    "\n",
    "fig, (axB,axA,axphase) = plt.subplots(1,3)\n",
    "\n",
    "axB.plot(sol1.t, sol1.y[0])\n",
    "axB.set_title('B ')\n",
    "axB.set(xlabel=\"t\")\n",
    "\n",
    "axA.plot(sol1.t, sol1.y[1])\n",
    "axA.set_title('A ')\n",
    "axA.set(xlabel=\"t\")\n",
    "\n",
    "axphase.plot(sol1.y[1], sol1.y[0])\n",
    "axphase.set_title('phaseplot')\n",
    "axphase.set(xlabel=\"A\", ylabel=\"B\")\n",
    "axphase.yaxis.set_label_position(\"right\")\n",
    "axphase.yaxis.tick_right()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bc875db5",
   "metadata": {},
   "source": [
    "## Solve system 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0dfb3a22",
   "metadata": {},
   "outputs": [],
   "source": [
    "# --------------------------------------------------------------------\n",
    "# Tragedy of the commons ODE system - course 25328 version\n",
    "# --------------------------------------------------------------------\n",
    "# state vector y = (B,A)\n",
    "# B: fish / km2\n",
    "# A: boats / km2\n",
    "# --- biological parameters ---\n",
    "# r: growth rate            [1/year]     \n",
    "# K: carrying capacity      [fish/km2]\n",
    "# b: habitat clearance rate [km2/boat/year] \n",
    "# --- economical parameters ---\n",
    "# s: sector adaptation rate [boats/km2/year]\n",
    "# p: catch price            [$/fish]       \n",
    "# c: operating cost         [$/year/boat]\n",
    "#\n",
    "# In this version parameter values reflect a semireal situation,\n",
    "# fishermen exploiting a fish population, and parameters/variables\n",
    "# are with units and strictly interpretable/observable, making them\n",
    "# easier to guesstimate. This is the version we use in course 25328 exercises.\n",
    "# It has 2 parameter settings - one with an\n",
    "# interior equilibrium, and one with a boundary equilibrium (sector crash)\n",
    "# By non-dimensionalization the system can be reduced to 3 essential\n",
    "# parameters. \n",
    "# --------------------------------------------------------------------\n",
    "def ToC(t,y, r,K,p,b,c,s):\n",
    "    B,A = y\n",
    "    dBdt = r*B*(1-B/K) - b*A*B\n",
    "    dAdt = s*(p*b*B/c - 1)\n",
    "    if A<0:\n",
    "        dAdt = max(dAdt, 0)     # dAdt > 0 is OK if A<0, allow fleet to rebuild\n",
    "    return np.array([dBdt,dAdt])\n",
    "\n",
    "# ---- this parameter combination corresponds to an interior equilibrium ----\n",
    "r  = 1            # [1/year]          growth rate\n",
    "K  = 1e3          # [fish/km2]        fish carrying capacity\n",
    "p  = 1            # [$/fish]          catch price\n",
    "b  = 50*5*1e-2*86 # [km2/boat/year]   habitat clearance rate (sail_days * fish_hours_pr_day * gear_width * sail_speed)   \n",
    "c  = 1e5          # [$/year/boat]     operating cost\n",
    "s  = 0.2/3/(3-1)  # [boats/km2/year]  sector adaptation rate (boat_dens / adapt_time / incentive_factor)\n",
    "y0 = (K, 0) # initial condition\n",
    "\n",
    "# ---- this parameter combination is in sector crash regime; uncomment line below\n",
    "# c   = 3e5; y0  = (1.3*K, 0.001) \n",
    "\n",
    "tspan  = (0, 30)\n",
    "t_eval = np.linspace(tspan[0], tspan[1], 1000)\n",
    "args   = (r,K,p,b,c,s)\n",
    "sol1   = solve_ivp(ToC, tspan, y0, t_eval=t_eval, args=args)\n",
    "\n",
    "# -------- plotting --------\n",
    "\n",
    "fig, (axB,axA,axphase) = plt.subplots(1,3)\n",
    "\n",
    "axB.plot(sol1.t, sol1.y[0])\n",
    "axB.set_title('B [fish/km2]')\n",
    "axB.set(xlabel=\"t [years]\")\n",
    "\n",
    "axA.plot(sol1.t, sol1.y[1])\n",
    "axA.set_title('A [boats/km2]')\n",
    "axA.set(xlabel=\"t [years]\")\n",
    "\n",
    "axphase.plot(sol1.y[1], sol1.y[0])\n",
    "axphase.set_title('phaseplot')\n",
    "axphase.set(xlabel=\"A [boats/km2]\", ylabel=\"B [fish/km2]\")\n",
    "axphase.yaxis.set_label_position(\"right\")\n",
    "axphase.yaxis.tick_right()\n",
    "\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "jupytext": {
   "text_representation": {
    "extension": ".md",
    "format_name": "myst",
    "format_version": 0.13,
    "jupytext_version": "1.17.2"
   }
  },
  "kernelspec": {
   "display_name": "InterMat",
   "language": "python",
   "name": "python3"
  },
  "source_map": [
   12,
   16,
   20,
   27,
   31,
   97,
   101
  ]
 },
 "nbformat": 4,
 "nbformat_minor": 5
}