{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Calculate probabilistic impact yearset\n",
    "\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_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 (given that the original impact object contains enough events). When setting `with_replacement=True`, the impact events are sampled with replacement . Note that sampling without replacement can lead to distorted sampling if the frequencies of the different impacts (`imp.frequency`) are not equal.\n",
    "\n",
    "*Correction factor.* By 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. Applying the correction factor can be avoided by setting `correction_fac=False`.\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. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "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([0, 1, 1, 0, 2, 3, 3, 0, 2, 3])"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 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(n_sampled_years, lam)\n",
    "events_per_year"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[array([], dtype=int64),\n",
       " array([0]),\n",
       " array([1]),\n",
       " array([], dtype=int64),\n",
       " array([4, 2]),\n",
       " array([2, 0, 4]),\n",
       " array([5, 1, 3]),\n",
       " array([], dtype=int64),\n",
       " array([5, 3]),\n",
       " array([5, 3, 1])]"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# generate the sampling vector\n",
    "sampling_vect = yearsets.sample_events(events_per_year, imp.frequency)\n",
    "sampling_vect"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0,  0,  2,  0, 64, 64, 70,  0, 68, 70])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# calculate the impact per year\n",
    "imp_per_year = yearsets.compute_imp_per_year(imp, sampling_vect)\n",
    "imp_per_year"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.7928994082840237"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# calculate the correction factor\n",
    "correction_factor = yearsets.calculate_correction_fac(imp_per_year, imp)\n",
    "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": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "yimp.at_event =  [  2   0   0  70 130   0   0 128  68  62]\n",
      "imp_per_year =  [  2   0   0  70 130   0   0 128  68  62]\n",
      "The expected annual impact 46.0 differs from the one of the original impact (26.8).\n"
     ]
    }
   ],
   "source": [
    "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, 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(\"yimp.at_event = \", yimp.at_event)\n",
    "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": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2025-09-19 15:39:54,743 - climada.util.yearsets - INFO - The correction factor is 0.5826086956521739.\n",
      "yimp.at_event =  [ 1.17  0.    0.   40.78 75.74  0.    0.   74.57 39.62 36.12]\n",
      "imp_per_year =  [ 1.17  0.    0.   40.78 75.74  0.    0.   74.57 39.62 36.12]\n",
      "The expected annual impact 26.8 is equal to the one of the original impact (26.8).\n"
     ]
    }
   ],
   "source": [
    "# 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 from imp_per_year computed above without 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",
    "    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": "climada_env600",
   "language": "python",
   "name": "python3"
  },
  "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.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}