update the tutorial

Author: ValentinGebhart

doc/user-guide/climada_util_yearsets.ipynb

@@ -6,13 +6,21 @@
    "source": [
     "# Calculate probabilistic impact yearset\n",
     "\n",
-    "This module generates a yearly impact `yimp` object which contains probabilistic annual impacts for a specified amount of years (`sampled_years`). The impact values are extracted from a given impact `imp` object that contains impact values per event. The amount of `sampled_years` can be specified as an integer or as a list of years to be sampled for. The amount of events per sampled year (`events_per_year`) are determined with a Poisson distribution centered around n_events per year (`lam` = sum(event_impacts.frequency). Then, the probabilistic events occurring in each sampled year are sampled uniformly from the input `imp` object and summed up per year. Thus, the `yimp` object contains the sum of sampled (event) impacts for each sampled year. In contrast to the expected annual impact (eai), an `yimp` object contains an impact for EACH sampled year and this value differs among years. The number of events_per_year and the selected_events are saved in a sampling vector (`sampling_vect`). \n",
+    "This module generates a yearly impact object `yimp` which contains probabilistic annual impacts for a specified amount of years (`sampled_years`). The impact values are extracted from a given impact `imp` object that contains impact values per event. The amount of sampled years as a list of years (`sampled_years`) to be sampled for. The amount of events per sampled year (`events_per_year`) are determined with a Poisson distribution (`lam` = sum(event_impacts.frequency)). Then, the events occurring in each sampled year are sampled from the input `imp` object and summed up per year. Thus, the `yimp` object contains the sum of sampled (event) impacts for each sampled year. In contrast to the expected annual impact (eai), an `yimp` object contains an impact for EACH sampled year and this value differs among years. The number of events_per_year and the selected_events are saved in a sampling vector (`sampling_vect`). \n",
     "\n",
-    "The function impact_yearsets performs all these computational steps, taking an `imp` and the number of sampled_years (`sampled_years`) as input. The output of the function is the `yimp` object and the `sampling_vect`.\n",
-    "Moreover, a `sampling_vect` (generated in a previous run) can be provided as optional input and the user can define `lam` and decide whether a correction factor shall be applied (the default is applying the correction factor).\n",
-    "Reapplying the same sampling_vect does not only allow to reproduce the generated `yimp`, but also for a physically consistent way of sampling impacts caused by different hazards. \n",
-    "The correction factor that is applied when the optional input `correction_fac`= True is a scaling of the computed `yimp` that assures that the eai(`yimp`) = eai(`imp`).\n",
+    "The function `impact_yearset` performs all these computational steps, taking an `imp` and the list of sampled_years (`sampled_years`) as input. The output of the function is the `yimp` object and the `sampling_vect`.\n",
+    "Moreover, a `sampling_vect` (generated in a previous run) can be provided as optional input and the user can custom-define the Poisson parameter `lam`. Reapplying the same sampling_vect does not only allow to reproduce the generated `yimp`, but also for a physically consistent way of sampling impacts caused by different hazards. \n",
     "\n",
+    "*Sampling options.* Per default, impact events are sampled without replacement (if the original impact object contains enough events), such that the different yearly impacts stem from different events. This is only implemented if the event frequencies of the original impact object (`imp.frequency`) are constant. If the events differ in their frequency, they are sampled with replacement. This sampling behaviour can also be achieved by setting `with_replacement=True`. \n",
+    "\n",
+    "*Correction factor.* Per default, a correction factor is applied uniformly to all yearly impacts, such that the final `yimp` object has the same average annual impact than `imp`, the original impact object. When setting `correction_fac=False`, the correction factor is not applied.\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
     "To make the process more transparent, this tutorial shows the single computations that are performed when generating an `yimp` object for a dummy event_impacts object. "
    ]
   },
@@ -22,60 +30,79 @@
    "metadata": {
     "scrolled": true
    },
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "import climada.util.yearsets as yearsets\n",
+    "from climada.engine import Impact\n",
+    "\n",
+    "# dummy event_impacts object containing 12 event impacts\n",
+    "imp = Impact(\n",
+    "    event_id=np.arange(6) + 10,\n",
+    "    event_name=np.arange(6) + 10,\n",
+    "    date=np.arange(6),\n",
+    "    coord_exp=np.array([[1, 2], [1.5, 2.5]]),\n",
+    "    eai_exp=np.array([13.4, 13.4]),\n",
+    "    at_event=np.array([0, 2, 4, 6, 60, 62]),\n",
+    "    frequency=np.full(6, 0.2),\n",
+    "    aai_agg=26.8,\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Step-by-step calculation"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
    "outputs": [
     {
      "data": {
       "text/plain": [
-       "array([2, 2, 2, 0, 4, 5, 4, 2, 3, 1])"
+       "array([1, 3, 0, 0, 2, 0, 1, 2, 1, 3])"
       ]
      },
-     "execution_count": 1,
+     "execution_count": 2,
      "metadata": {},
      "output_type": "execute_result"
     }
    ],
    "source": [
-    "import numpy as np\n",
-    "\n",
-    "import climada.util.yearsets as yearsets\n",
-    "from climada.engine import Impact\n",
-    "\n",
-    "# dummy event_impacts object containing 10 event_impacts with the values 10-110\n",
-    "# and the frequency 0.2 (Return period of 5 years)\n",
-    "imp = Impact()\n",
-    "imp.at_event = np.arange(10, 110, 10)\n",
-    "imp.frequency = np.array(np.ones(10) * 0.2)\n",
-    "\n",
-    "# the number of years to sample impacts for (length(yimp.at_event) = sampled_years)\n",
-    "sampled_years = 10\n",
+    "# the number of years to sample impacts for (length(yimp.at_event) = n_sampled_years)\n",
+    "n_sampled_years = 10\n",
     "\n",
     "# sample number of events per sampled year\n",
     "lam = np.sum(imp.frequency)\n",
-    "events_per_year = yearsets.sample_from_poisson(sampled_years, lam)\n",
+    "events_per_year = yearsets.sample_from_poisson(n_sampled_years, lam)\n",
     "events_per_year"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 2,
+   "execution_count": 3,
    "metadata": {},
    "outputs": [
     {
      "data": {
       "text/plain": [
-       "[array([8, 3]),\n",
-       " array([7, 0]),\n",
-       " array([4, 6]),\n",
-       " array([], dtype=int32),\n",
-       " array([5, 9, 1, 2]),\n",
-       " array([1, 6, 0, 7, 2]),\n",
-       " array([4, 9, 5, 8]),\n",
-       " array([9, 8]),\n",
-       " array([5, 3, 4]),\n",
-       " array([1])]"
+       "[array([2]),\n",
+       " array([0, 4, 3]),\n",
+       " array([], dtype=int64),\n",
+       " array([], dtype=int64),\n",
+       " array([5, 1]),\n",
+       " array([], dtype=int64),\n",
+       " array([3]),\n",
+       " array([5, 1]),\n",
+       " array([0]),\n",
+       " array([2, 1, 5])]"
       ]
      },
-     "execution_count": 2,
+     "execution_count": 3,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -88,16 +115,16 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 3,
+   "execution_count": 4,
    "metadata": {},
    "outputs": [
     {
      "data": {
       "text/plain": [
-       "[130, 90, 120, 0, 210, 210, 300, 190, 150, 20]"
+       "array([ 4, 66,  0,  0, 64,  0,  6, 64,  0, 68])"
       ]
      },
-     "execution_count": 3,
+     "execution_count": 4,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -110,16 +137,16 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 4,
+   "execution_count": 5,
    "metadata": {},
    "outputs": [
     {
      "data": {
       "text/plain": [
-       "0.7746478873239436"
+       "0.9852941176470589"
       ]
      },
-     "execution_count": 4,
+     "execution_count": 5,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -130,65 +157,111 @@
     "correction_factor"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Direct calculation"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "If we use the same seed, the `yimp.at_event` values coincide with the step-by-step `imp_per_year` values, if no correction factor is applied."
+   ]
+  },
   {
    "cell_type": "code",
-   "execution_count": 5,
+   "execution_count": 6,
    "metadata": {},
    "outputs": [
     {
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "The yimp.at_event values equal our step-by-step computed imp_per_year:\n",
-      "yimp.at_event =  [90, 240, 150, 70, 40, 90, 60, 90, 170, 110]\n",
-      "imp_per_year =  [130, 90, 120, 0, 210, 210, 300, 190, 150, 20]\n"
+      "yimp.at_event =  [ 60   0   0   8 128   0   0 130 122   8]\n",
+      "imp_per_year =  [ 60   0   0   8 128   0   0 130 122   8]\n",
+      "The expected annual impact 45.6 differs from the one of the original impact (26.8).\n"
      ]
     }
    ],
    "source": [
-    "# compare the resulting yimp with our step-by-step computation without applying the correction factor:\n",
+    "example_seed = 12345  # fix a seed\n",
+    "sampled_years = list(range(2015, 2025))\n",
+    "\n",
+    "# direct computation\n",
     "yimp, sampling_vect = yearsets.impact_yearset(\n",
-    "    imp, sampled_years=list(range(1, 11)), correction_fac=False\n",
+    "    imp, sampled_years, correction_fac=False, seed=example_seed\n",
+    ")\n",
+    "\n",
+    "# step-by-step computation without applying the correction factor\n",
+    "events_per_year = yearsets.sample_from_poisson(\n",
+    "    len(sampled_years), lam, seed=example_seed\n",
+    ")\n",
+    "sampling_vect = yearsets.sample_events(\n",
+    "    events_per_year, imp.frequency, seed=example_seed\n",
     ")\n",
+    "imp_per_year = yearsets.compute_imp_per_year(imp, sampling_vect)\n",
+    "\n",
     "\n",
-    "print(\"The yimp.at_event values equal our step-by-step computed imp_per_year:\")\n",
     "print(\"yimp.at_event = \", yimp.at_event)\n",
-    "print(\"imp_per_year = \", imp_per_year)"
+    "print(\"imp_per_year = \", imp_per_year)\n",
+    "print(\n",
+    "    f\"The expected annual impact {round(yimp.aai_agg,2)} differs from the one of the original impact ({round(imp.aai_agg,2)}).\"\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Similarly, if we use the same seed, the `yimp.at_event` values when using the correction factor coincide with the step-by-step `imp_per_year` values after applyng the correction factor."
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 6,
+   "execution_count": 7,
    "metadata": {},
    "outputs": [
     {
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "The same can be shown for the case of applying the correction factor.The yimp.at_event values equal our step-by-step computed imp_per year:\n",
-      "yimp.at_event =  [ 54.54545455  47.72727273  95.45454545 109.09090909   0.\n",
-      "  40.90909091  27.27272727   0.         109.09090909  27.27272727]\n",
-      "imp_per_year =  [167.81818182 116.18181818 154.90909091   0.         271.09090909\n",
-      " 271.09090909 387.27272727 245.27272727 193.63636364  25.81818182]\n"
+      "2025-07-17 14:56:28,995 - climada.util.yearsets - INFO - The correction factor is 0.5877192982456141.\n",
+      "yimp.at_event =  [35.26  0.    0.    4.7  75.23  0.    0.   76.4  71.7   4.7 ]\n",
+      "imp_per_year =  [35.26  0.    0.    4.7  75.23  0.    0.   76.4  71.7   4.7 ]\n",
+      "The expected annual impact 26.8 is equal to the one of the original impact (26.8).\n"
      ]
     }
    ],
    "source": [
-    "# and here the same comparison with applying the correction factor (default settings):\n",
-    "yimp, sampling_vect = yearsets.impact_yearset(imp, sampled_years=list(range(1, 11)))\n",
+    "# direct computation\n",
+    "yimp, sampling_vect = yearsets.impact_yearset(\n",
+    "    imp, sampled_years, correction_fac=True, seed=example_seed\n",
+    ")\n",
     "\n",
+    "# compute correction factor\n",
+    "correction_factor = yearsets.calculate_correction_fac(imp_per_year, imp)\n",
+    "\n",
+    "print(\"yimp.at_event = \", np.round(yimp.at_event, 2))\n",
+    "print(\"imp_per_year = \", np.round(imp_per_year * correction_factor, 2))\n",
     "print(\n",
-    "    \"The same can be shown for the case of applying the correction factor.\"\n",
-    "    \"The yimp.at_event values equal our step-by-step computed imp_per year:\"\n",
-    ")\n",
-    "print(\"yimp.at_event = \", yimp.at_event)\n",
-    "print(\"imp_per_year = \", imp_per_year / correction_factor)"
+    "    f\"The expected annual impact {round(yimp.aai_agg,2)} is equal to the one of the original impact ({round(imp.aai_agg,2)}).\"\n",
+    ")"
    ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
   }
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python 3 (ipykernel)",
+   "display_name": "climada_env600",
    "language": "python",
    "name": "python3"
   },
@@ -202,7 +275,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.8.12"
+   "version": "3.11.11"
   }
  },
  "nbformat": 4,