{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "34a69e70-7601-4399-80a7-750e68829666",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "from tqdm import tqdm\n",
    "import egttools as egt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5ac97f32-5bbe-46e1-8317-f4f2a536f684",
   "metadata": {},
   "outputs": [],
   "source": [
    "nb_strategies = 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ecf8e325-8b7c-4dc3-861e-b7a209891827",
   "metadata": {},
   "outputs": [],
   "source": [
    "markers = [\"P\", \"o\", \"v\", \"h\", \"s\", \"X\", \"D\"]\n",
    "colors = sns.color_palette(\"icefire\", nb_strategies)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "83a3b29e-ed3a-4eef-9edb-bb2dfb929d3f",
   "metadata": {},
   "source": [
    "# Exercise 1: Replicator Dynamics of Sin City"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "de26f0a4-a15f-43f8-a281-836e0a5457a5",
   "metadata": {},
   "source": [
    "## Part 1: Analytical/Numerical calculations"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "22f98ae8-624a-45ab-a123-4dd0bec1ae91",
   "metadata": {},
   "source": [
    "### 1.1.1"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3c29204d-1d22-4edd-9907-53cbceb7c579",
   "metadata": {},
   "source": [
    "Write here your solution"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5ea961bb-3434-4ed0-9fbb-aae119490288",
   "metadata": {},
   "source": [
    "### 1.1.2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "34f1d7fc-4cc5-4d03-9932-b5477644ef81",
   "metadata": {},
   "outputs": [],
   "source": [
    "# define the parameters\n",
    "p_o = 0.2\n",
    "p_c = 0.5\n",
    "e_c = 1\n",
    "e_o = 2\n",
    "b_c = 15\n",
    "b_cf = 2\n",
    "b_o = 10\n",
    "b_p = 10\n",
    "c_cf = 5\n",
    "c_c = 50\n",
    "c_o = 3\n",
    "c_po = 1\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8fe9afcd-2255-4b44-83f0-89cf7e2f76a8",
   "metadata": {},
   "outputs": [],
   "source": [
    "strategy_labels = [\"O\", \"C\", \"P\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bbc76a0c-75aa-441d-9fb0-e774645f8b7c",
   "metadata": {},
   "source": [
    "## 1.1.3"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1b2aad2c-6eef-4f55-8dd3-b2b325d6390d",
   "metadata": {},
   "source": [
    "Before you start the next part, please assign the payoff matrix you have created in this part to a variable called `payoff_matrix_1`. This will allow you to run the example code we provide."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8ccb5fc4-80b2-4a5e-88dc-165fffbe0916",
   "metadata": {},
   "outputs": [],
   "source": [
    "payoff_matrix_1 = # assign a payoff matrix to this variable !"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b1663b0a-aefa-4846-9807-6c255981f88f",
   "metadata": {},
   "source": [
    "## Part 2 Monte Carlo simulations"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f871eb6f-3d1a-4451-be21-ddb014abeb01",
   "metadata": {},
   "source": [
    "First let's plot the time evolution of the simulation by integrating the differential equation with odeint in the time interval [0, 10] in 1000 evenly distributed steps over 1000 independent runs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2c712915-b17e-4a8b-a3e4-d1b4b40e62f4",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.integrate import odeint"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "17e91d6b-8b7f-430f-8df5-65fa0a62a727",
   "metadata": {},
   "outputs": [],
   "source": [
    "nb_runs = 1000\n",
    "nb_time_steps = 1000\n",
    "t = np.arange(0, 10, 10/nb_time_steps)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7519b173-d292-411f-89a7-ac17640b9054",
   "metadata": {},
   "outputs": [],
   "source": [
    "def run_replicator_dynamics(t: np.ndarray, payoff_matrix: np.ndarray, nb_runs: int):\n",
    "    results = []\n",
    "    for i in range(nb_runs):\n",
    "        x0 = egt.sample_unit_simplex(3)\n",
    "        result = odeint(lambda x, t: egt.analytical.replicator_equation(x, payoff_matrix), x0, t)\n",
    "        results.append(result)\n",
    "        \n",
    "    return results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6ac9c3b6-b29d-4ce3-81ee-f6a1d6bb6ae5",
   "metadata": {},
   "outputs": [],
   "source": [
    "results = run_replicator_dynamics(t, payoff_matrix_1, nb_runs)\n",
    "\n",
    "results = np.asarray(results)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6a3fcdfe-1693-457f-a3a3-c8c47799db96",
   "metadata": {},
   "source": [
    "The following plot shows that at $t=4$ most of the simulations have already converged. This is important, since it means that all simulations converge before $t=10$ and thus our avarages the following exercise will be correct. The convergence, as expected, depends on the starting point, and on average, the population will end more often with 100% of Poice Officers (O). This means that, as could be seen in the simplex plot you should produce for Part 1, O has the largest basin of attraction in this configuration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "086b95d6-3d2b-40a1-92b0-f634ef731b98",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(10, 3))\n",
    "\n",
    "for run in results:\n",
    "    for i in range(nb_strategies):\n",
    "        ax.plot(t, run[:, i],\n",
    "            linewidth=.05, alpha=0.6, color=colors[i])\n",
    "        \n",
    "for i in range(nb_strategies):\n",
    "    ax.plot(t, np.mean(results[:, :, i], axis=0), linewidth=1.5, \n",
    "            alpha=1, color=colors[i], label=strategy_labels[i])   \n",
    "\n",
    "ax.legend(frameon=False, bbox_to_anchor=(1.1,1,0,0), loc='upper right')\n",
    "ax.set_ylabel(\"frequency\", fontsize=14)\n",
    "ax.set_xlabel(\"time step, $t$\", fontsize=14)\n",
    "# ax.set_ylim(-0.2, 1.2)\n",
    "sns.despine()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dd5aff65-316f-4176-a7b3-b4efdff2dc2d",
   "metadata": {},
   "source": [
    "### 1.2.1 and 1.2.2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "455303fb-9175-456f-9162-1a9de6dda90a",
   "metadata": {},
   "source": [
    "Now let's apply the same methodology to run the Monte-Carlo simulations required to solve these two exercises"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "703b534e-b161-4ca3-b47d-baf4bfa904bc",
   "metadata": {},
   "outputs": [],
   "source": [
    "p_os = np.linspace(0, 1, 11)\n",
    "nb_runs = 1000\n",
    "nb_generations = 1000\n",
    "t = np.arange(0, 10, 10/nb_generations)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "85105075-6753-4491-bdbc-38e4fb8a8b69",
   "metadata": {},
   "source": [
    "You can use the following code to produce the plot with the results of this exercise (this is optional). Adapt the parameter names to whatever you are using."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f98ad3fb-5de1-471d-ba2c-83765200e7d8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Code to produce the plot \n",
    "fig, ax = plt.subplots(figsize=(10, 3))\n",
    "\n",
    "for i in range(nb_strategies):\n",
    "    ax.plot(p_os, avg_results[:, i], marker=markers[i], \n",
    "            label=strategy_labels[i], color=colors[i], lw=2)\n",
    "\n",
    "ax.legend(frameon=False)\n",
    "ax.set_ylabel(\"frequency\", fontsize=14)\n",
    "ax.set_xlabel(\"probability of catching a criminal, $p_o$\", fontsize=14)\n",
    "# ax.set_ylim(-0.2, 1.2)\n",
    "sns.despine()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4eb9f7b1-ea91-4b7f-9c6c-bb4137db9610",
   "metadata": {},
   "source": [
    "### 1.2.3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e93d3d7e-3e52-4fc3-8f4c-15d8acb51ee2",
   "metadata": {},
   "outputs": [],
   "source": [
    "p_cs = np.linspace(0, 1, 11)\n",
    "nb_runs = 1000\n",
    "nb_generations = 1000\n",
    "p_o = 0.2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c8897f14-1be9-4643-8b58-192ac1fabf81",
   "metadata": {},
   "source": [
    "You can use the following code to produce the plot with the results of this exercise (this is optional). Adapt the parameter names to whatever you are using."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "392164ef-2404-4985-8c2d-834db5c91eb8",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(10, 3))\n",
    "\n",
    "for i in range(nb_strategies):\n",
    "    ax.plot(p_os, avg_results[:, i], marker=markers[i], \n",
    "            label=strategy_labels[i], color=colors[i], lw=2)\n",
    "    \n",
    "ax.legend(frameon=False)\n",
    "ax.set_ylabel(\"frequency\", fontsize=14)\n",
    "ax.set_xlabel(\"probability a criminal effectively robs a common citizen, $p_c$\", fontsize=14)\n",
    "# ax.set_ylim(-0.2, 1.2)\n",
    "sns.despine()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a56e37e4-3c3a-487c-8a72-f3bfad7bd45c",
   "metadata": {},
   "source": [
    "# Exercise 2: Stochastic Dynamics of Sin City"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6bc384c4-52e0-445f-b57b-3ae393ee0df3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# reinitialize the changed values\n",
    "p_o = 0.2\n",
    "p_c = 0.5"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c3f894a1-5269-4572-bbc9-d260962942b5",
   "metadata": {},
   "source": [
    "First, as we did in Part 2 of Exercise 1, let's plot the time evolution of the simulation over 1000 time-steps in 1000 independent runs. However, since now we have a stochastic model, we can no longer use ode-int to integrate. Instead, let's use the PairwiseComparisonNumerical class of `egttools` to evolve the population. Since we now want to see what happens in each time step, we will use the method `run`. However, for the next exercises, all you are interested in is the last state the population reaches, so you should use the method `evolve` instead. The call to `evolve` requires the following parameters `evolve(nb_generations, beta, mu, init_state)` where `nb_generations` is the number of generations for which the simulation will be run, `beta` and `mu` are the parameters of the model, and `init_state` is the initial state of the population."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "80129247-1cfa-43e2-9b3e-58c58b6df7bc",
   "metadata": {},
   "source": [
    "We will first instantiate a game class called `Matrix2PlayerGameHolder` that will allow us to pass our payoff matrix to the `PairwiseComparisonNumerical` class, which will run the simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "193411d9-0af3-4299-9f86-b293559a82e0",
   "metadata": {},
   "outputs": [],
   "source": [
    "from egttools.games import Matrix2PlayerGameHolder\n",
    "from egttools.numerical import PairwiseComparisonNumerical"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "94d64768-35bf-4c05-8d50-90225fb83045",
   "metadata": {},
   "source": [
    "Then we define the parameters that are related to the Stochastic model of Social Learning (cache_size is a parameter of the class PairwiseComparisonNumerical that defines how much cache memory will be available, this accelerates computations)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f27b5212-f24f-424e-aa5e-59a684de4129",
   "metadata": {},
   "outputs": [],
   "source": [
    "Z = 100\n",
    "beta = 1\n",
    "mu = 1e-3\n",
    "cache_size = 100000\n",
    "nb_population_states = egt.calculate_nb_states(Z, nb_strategies)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "73c67143-de00-41fa-9cf4-a0843c465c3a",
   "metadata": {},
   "source": [
    "Then we define how many runs and generations we will compute. We will continue to use in this example the payoff matrix computed in the first exercise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d1c6a2ad-18ae-4436-8280-83403c903b6c",
   "metadata": {},
   "outputs": [],
   "source": [
    "nb_runs = 1000\n",
    "nb_generations = 1000"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "73673f38-5394-4a70-9114-4a1d9c25475a",
   "metadata": {},
   "source": [
    "Now we instationate the game and the evolver. Bear in mind for the next exercise that you will need to instantiate a new game and evolver for every new payoff matrix you calculate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "527c8643-6dac-47c8-aecb-19c3c48632c4",
   "metadata": {},
   "outputs": [],
   "source": [
    "game = Matrix2PlayerGameHolder(nb_strategies, payoff_matrix_1)\n",
    "evolver = PairwiseComparisonNumerical(Z, game, cache_size)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "edeccfd5-5684-43b6-ba26-bd1a9bb0d07f",
   "metadata": {},
   "outputs": [],
   "source": [
    "results = []\n",
    "\n",
    "for i in range(nb_runs):\n",
    "    index = np.random.randint(0, nb_population_states)\n",
    "    x0 = egt.sample_simplex(index, Z, nb_strategies)\n",
    "    result = evolver.run(nb_generations, 0, beta, mu, x0)\n",
    "    results.append(result)\n",
    "results = np.asarray(results)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e4d81abb-97cc-4669-b91e-554c0f43d306",
   "metadata": {},
   "source": [
    "There are a few differences in the following figure with respect of the model with infinite populations (which uses the replicator equation) which you should notice: \n",
    "\n",
    "1. the population now varies in discrete steps of at most 1 individual (since we have a birth-death process)\n",
    "2. instead of frequency the state of the population is now represented by the counts of each strategy in the finite population (or by its proportion).\n",
    "3. many simulations have not yet converged by generation 1000."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f2c131d9-1b63-4f51-8f1c-188255837cf5",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(10, 3))\n",
    "\n",
    "for run in results:\n",
    "    for i in range(nb_strategies):\n",
    "        ax.plot(range(nb_generations), run[:, i]/Z,\n",
    "            linewidth=.05, alpha=0.6, color=colors[i])\n",
    "        \n",
    "for i in range(nb_strategies):\n",
    "    ax.plot(range(nb_generations), np.mean(results[:, :, i]/Z, axis=0), linewidth=1.5, \n",
    "            alpha=1, color=colors[i], label=strategy_labels[i])   \n",
    "\n",
    "ax.legend(frameon=False, bbox_to_anchor=(1.1,1,0,0), loc='upper right')\n",
    "ax.set_ylabel(\"proportion ($k/Z$)\", fontsize=14)\n",
    "ax.set_xlabel(\"generation\", fontsize=14)\n",
    "# ax.set_ylim(-0.2, 1.2)\n",
    "sns.despine()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6f7979f2-8647-4868-8aeb-46426cae2494",
   "metadata": {},
   "source": [
    "Point 3 of the previous discussion above the previous plot indicates that we need to run more generations of the simulation to have reliable results. Let's do that:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4276a5d0-25e2-431b-a5e8-b356cedcf735",
   "metadata": {},
   "outputs": [],
   "source": [
    "nb_runs = 1000\n",
    "nb_generations = 10000 # now we use 10ˆ4 generations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "36bf8128-31bf-46b1-85d1-8478c88b4fd1",
   "metadata": {},
   "outputs": [],
   "source": [
    "results = []\n",
    "\n",
    "for i in range(nb_runs):\n",
    "    index = np.random.randint(0, nb_population_states)\n",
    "    x0 = egt.sample_simplex(index, Z, nb_strategies)\n",
    "    result = evolver.run(nb_generations, 0, beta, mu, x0)\n",
    "    results.append(result)\n",
    "results = np.asarray(results)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c05edc49-43f8-4582-ba21-c015f11f3e0f",
   "metadata": {},
   "source": [
    "The following plot shows that $10^4$ generations should be enough, and this is the number of generations that you need to use for Exercise 2. However, you should notice that there are some ciclic dynamics occuring until very late in the simulations, and there are small variations even at the end. If there was no mutation, we would not observe those variations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1c6141c8-7e5c-436c-8f06-2c62424224b1",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(10, 3))\n",
    "\n",
    "for run in results:\n",
    "    for i in range(nb_strategies):\n",
    "        ax.plot(range(nb_generations), run[:, i]/Z,\n",
    "            linewidth=.05, alpha=0.6, color=colors[i])\n",
    "        \n",
    "for i in range(nb_strategies):\n",
    "    ax.plot(range(nb_generations), np.mean(results[:, :, i]/Z, axis=0), linewidth=1.5, \n",
    "            alpha=1, color=colors[i], label=strategy_labels[i])   \n",
    "\n",
    "ax.legend(frameon=False, bbox_to_anchor=(1.1,1,0,0), loc='upper right')\n",
    "ax.set_ylabel(\"proportion ($k/Z$)\", fontsize=14)\n",
    "ax.set_xlabel(\"generation\", fontsize=14)\n",
    "# ax.set_ylim(-0.2, 1.2)\n",
    "sns.despine()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a174b61e-43eb-4a1d-b071-5d185a96d4b1",
   "metadata": {},
   "source": [
    "### 2.1.1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3507fbcf-7241-400e-b806-5aa506ba7903",
   "metadata": {},
   "outputs": [],
   "source": [
    "p_os = np.linspace(0, 1, 11)\n",
    "nb_runs = 1000\n",
    "nb_generations = 10000"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4acb2b27-3e71-4f3e-8a40-5c64502324aa",
   "metadata": {},
   "source": [
    "You can use the following code to produce the plot with the results of this exercise (this is optional). Adapt the parameter names to whatever you are using."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3212e184-f7c6-4644-812a-0af97de82b31",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(10, 3))\n",
    "\n",
    "for i in range(nb_strategies):\n",
    "    ax.plot(p_os, avg_results[:, i], marker=markers[i], \n",
    "            label=strategy_labels[i], color=colors[i], lw=2)\n",
    "\n",
    "ax.legend(frameon=False)\n",
    "ax.set_ylabel(\"frequency\", fontsize=14)\n",
    "ax.set_xlabel(\"probability of catching a criminal, $p_o$\", fontsize=14)\n",
    "# ax.set_ylim(-0.2, 1.2)\n",
    "sns.despine()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f937df3d-a075-4e74-b2ce-7a66a93c56cb",
   "metadata": {
    "tags": []
   },
   "source": [
    "### 2.1.2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "023c7b57-96c4-4f67-b020-3181fadb545d",
   "metadata": {},
   "outputs": [],
   "source": [
    "p_cs = np.linspace(0, 1, 11)\n",
    "nb_runs = 1000\n",
    "nb_generations = 10000\n",
    "p_o = 0.2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6ee56c35-b7e5-4ebf-9553-8786bf77501c",
   "metadata": {},
   "source": [
    "You can use the following code to produce the plot with the results of this exercise (this is optional). Adapt the parameter names to whatever you are using."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b57abab2-23ac-428a-9b76-05ad3c23a23a",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(10, 3))\n",
    "\n",
    "for i in range(nb_strategies):\n",
    "    ax.plot(p_os, avg_results[:, i], marker=markers[i], \n",
    "            label=strategy_labels[i], color=colors[i], lw=2)\n",
    "    \n",
    "ax.legend(frameon=False)\n",
    "ax.set_ylabel(\"proportion ($k/Z$)\", fontsize=14)\n",
    "ax.set_xlabel(\"probability a criminal effectively robs a common citizen, $p_c$\", fontsize=14)\n",
    "# ax.set_ylim(-0.2, 1.2)\n",
    "sns.despine()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c68336ea-ca5c-4f13-a04c-ccff25afcb4e",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "python3 (trust_dynamics)",
   "language": "python",
   "name": "conda-env-trust_dynamics-py"
  },
  "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.11.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}