{ "cells": [ { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "# Uncertainty assessment\n", "\n", "The IPCC guidelines provide information on taking into account uncertainty.\n", "The `bonsai_ipcc` package implements those by allowing the user to do analytical error propagation and Monte Carlo simualtion when running a tier sequence. \n", "\n", "To use one of the approaches, each parameter involved in a sequence requires information about its uncertainty. This information is specifeid by the properties \"def\" , \"min\", \"max\", \"abs_min\" and \"abs_max\" (mean, 2.5 percentile, 97.5 percentile, absolute minimum, absolute maximum) in the parameter table.\n", "\n", "Let´s use have a look into the waste incineration example." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Add data with uncertainty\n", "As already done in the previous tutorial for running a sequence, we need to add the data which is not provided by the package. So let´s do this again." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import bonsai_ipcc\n", "import pandas as pd\n", "import numpy as np\n", "\n", "my_ipcc = bonsai_ipcc.IPCC()\n", "\n", "# urban population\n", "d = {\n", " \"year\": [2010,2010,2010,2010,2010],\n", " \"region\": [\"DE\",\"DE\",\"DE\",\"DE\",\"DE\"],\n", " \"property\": [\n", " \"def\",\"min\",\"max\",\"abs_min\",\"abs_max\"\n", " ],\n", " \"value\": [\n", " 62940432.0,61996325.52,63884538.48,0.0,np.inf,\n", " ],\n", " \"unit\": [\n", " \"cap/yr\",\"cap/yr\",\"cap/yr\",\"cap/yr\",\"cap/yr\",\n", " ],\n", "}\n", "urb_pop = pd.DataFrame(d).set_index([\"year\", \"region\", \"property\"])\n", "\n", "my_ipcc.waste.incineration.parameter.urb_population=urb_pop" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To take the uncertainty into account, we need to add the values for `def`, `min`, `max`, `abs_min`, and `abs_max`. A description of the `property` dimension can be found here:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
descriptionremarks
code
defdefaultmean
minminimum2.5th percentile
maxmaximum97.5th percentile
abs_maxabsolute maximuntheoretical upper bound
abs_minabsolute minimumtheoretical lower bound
dummyplace holder for emtpy tablesNaN
\n", "
" ], "text/plain": [ " description remarks\n", "code \n", "def default mean\n", "min minimum 2.5th percentile\n", "max maximum 97.5th percentile\n", "abs_max absolute maximun theoretical upper bound\n", "abs_min absolute minimum theoretical lower bound\n", "dummy place holder for emtpy tables NaN" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "my_ipcc.waste.incineration.dimension.property" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## Monte Carlo simulation\n", "\n", "To run the sequence as Monte Carlo simulation you can use `monte_carlo` as keyword. Notably, the ipcc package uses a vectorized way of iplmenting Monte Carlo simulation. Thus, all parameters for those uncertainty information is included in the parameter tables and all calculated results are represented as numpy arrays (1000 values).\n", "\n", "We use the incineration example again, but chose the paper fraction of municipal waste instead of plastics." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "2023-12-19 15:35:42,594 - INFO - Incineration sequence started --->\n", "2023-12-19 15:35:42,595 - INFO - Uncertainty distribution for parameter 'urb_population':\n", "2023-12-19 15:35:42,595 - INFO - lognormal distribution\n", "/Users/TN76JP/Documents/coderefinery/ipcc/ipcc/src/ipcc/_sequence.py:87: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`\n", " new_c = conc_df.loc[c][j]\n", "2023-12-19 15:35:42,598 - INFO - Uncertainty distribution for parameter 'msw_gen_rate':\n", "2023-12-19 15:35:42,598 - INFO - lognormal distribution\n", "2023-12-19 15:35:42,598 - INFO - 'Coordinates (2010, 'DE')' has been replaced by '[2006, 'Western Europe']' during reading parameter table 'msw_gen_rate'\n", "2023-12-19 15:35:42,599 - INFO - Uncertainty distribution for parameter 'msw_frac_to_incin':\n", "2023-12-19 15:35:42,600 - INFO - normal distribution, lower uncertainty\n", "2023-12-19 15:35:42,600 - INFO - 'Coordinates (2010, 'DE')' has been replaced by '[2006, 'Western Europe']' during reading parameter table 'msw_frac_to_incin'\n", "2023-12-19 15:35:42,602 - INFO - Uncertainty distribution for parameter 'msw_type_frac':\n", "2023-12-19 15:35:42,603 - INFO - normal distribution, lower uncertainty\n", "2023-12-19 15:35:42,603 - INFO - 'Coordinates (2010, 'DE', 'msw_paper')' has been replaced by '[2006, 'DE', 'msw_paper']' during reading parameter table 'msw_type_frac'\n", "2023-12-19 15:35:42,605 - INFO - Uncertainty distribution for parameter 'incintype_frac':\n", "/Users/TN76JP/Documents/coderefinery/ipcc/ipcc/src/ipcc/uncertainties.py:295: RuntimeWarning: invalid value encountered in double_scalars\n", " if max(sd / (mean - abs_min), sd / (abs_max - mean)) <= 0.3:\n", "2023-12-19 15:35:42,606 - INFO - normal distribution, lower uncertainty\n", "2023-12-19 15:35:42,606 - INFO - 'Coordinates (2010, 'DE', 'inc_unspecified')' has been replaced by '[2006, 'World', 'inc_unspecified']' during reading parameter table 'incintype_frac'\n", "2023-12-19 15:35:42,609 - INFO - Uncertainty distribution for parameter 'dm':\n", "2023-12-19 15:35:42,609 - INFO - truncated normal distribution with adjusting based on Rodriques 2015 (moderate)\n", "2023-12-19 15:35:42,611 - INFO - 'Coordinates (2010, 'DE', 'msw_paper')' has been replaced by '[2006, 'World', 'msw_paper']' during reading parameter table 'dm'\n", "2023-12-19 15:35:42,612 - INFO - Uncertainty distribution for parameter 'cf':\n", "2023-12-19 15:35:42,613 - INFO - normal distribution, lower uncertainty\n", "2023-12-19 15:35:42,613 - INFO - 'Coordinates (2010, 'DE', 'msw_paper')' has been replaced by '[2006, 'World', 'msw_paper']' during reading parameter table 'cf'\n", "2023-12-19 15:35:42,614 - INFO - Uncertainty distribution for parameter 'fcf':\n", "2023-12-19 15:35:42,614 - INFO - truncated normal distribution with adjusting based on Rodriques 2015 (moderate)\n", "2023-12-19 15:35:42,615 - INFO - 'Coordinates ('DE', 'msw_paper', 'inc_unspecified')' has been replaced by '['World', 'msw_paper', 'inc_unspecified']' during reading parameter table 'fcf'\n", "2023-12-19 15:35:42,616 - INFO - Uncertainty distribution for parameter 'of':\n", "2023-12-19 15:35:42,616 - INFO - normal distribution, lower uncertainty\n", "2023-12-19 15:35:42,617 - INFO - 'Coordinates ('DE', 'msw_paper', 'inc_unspecified')' has been replaced by '['World', 'msw_paper', 'inc_unspecified']' during reading parameter table 'of'\n", "2023-12-19 15:35:42,617 - INFO - ---> Incineration sequence finalized.\n" ] } ], "source": [ "\n", "\n", "# run the sequence\n", "my_tier = my_ipcc.waste.incineration.sequence.tier1_co2(year=2010, region=\"DE\", wastetype= \"msw_paper\", incintype= \"inc_unspecified\", uncertainty=\"monte_carlo\")\n", "s = my_tier.to_dict()\n" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "numpy.ndarray" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# check the last step (emissions)\n", "type(s[\"co2_emissions\"].value)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Plot the results" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(array([119., 229., 228., 183., 107., 63., 41., 20., 3., 7.]),\n", " array([1.02748936e-01, 2.78783377e+01, 5.56539265e+01, 8.34295153e+01,\n", " 1.11205104e+02, 1.38980693e+02, 1.66756282e+02, 1.94531871e+02,\n", " 2.22307459e+02, 2.50083048e+02, 2.77858637e+02]),\n", " )" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAAsTAAALEwEAmpwYAAAM6klEQVR4nO3dX4xc9XnG8e9TSLhIkIB6a1lguiTyjXtRYq0oUlFEhZSAuTC5QXBRrAjJvQApkdoLp7kIN5GcSklVpBTJKSimSqFICcIStA21IqFekMRExPwrwU2MsGVsp1SEKlJayNuLPW4my673z+x6PC/fjzSaM79zZs/76qyePfObmbOpKiRJvfzOpAuQJK0/w12SGjLcJakhw12SGjLcJamhiyddAMCmTZtqdnZ20mVI0lR57rnnfl5VM4utuyDCfXZ2lsOHD0+6DEmaKkleX2qd0zKS1JDhLkkNGe6S1JDhLkkNGe6S1JDhLkkNGe6S1JDhLkkNGe6S1NAF8Q3VaTW798mJ7PfYvlsnsl9J08Nwn0KT+qMC/mGRpoXTMpLUkOEuSQ0Z7pLUkOEuSQ0Z7pLUkOEuSQ0Z7pLUkOEuSQ0Z7pLUkOEuSQ0Z7pLUkOEuSQ0Z7pLUkOEuSQ0Z7pLUkOEuSQ0Z7pLUkOEuSQ0Z7pLUkOEuSQ0Z7pLU0LLhnmRrku8leTnJS0k+N4xfkeTpJK8N95cP40lyf5KjSY4k2bHRTUiSfttKztzfBf68qrYD1wP3JNkO7AUOVdU24NDwGOAWYNtw2wM8sO5VS5LOadlwr6qTVfWjYfkd4BXgSmAXcGDY7ABw27C8C3i45j0LXJZky3oXLkla2qrm3JPMAp8Avg9srqqTw6o3gc3D8pXAGyNPOz6MLfxZe5IcTnL4zJkzq61bknQOKw73JB8Fvg18vqp+Mbquqgqo1ey4qvZX1VxVzc3MzKzmqZKkZawo3JN8iPlg/1ZVfWcYPnV2umW4Pz2MnwC2jjz9qmFMknSerOTTMgEeBF6pqq+NrDoI7B6WdwNPjIzfNXxq5nrg7ZHpG0nSeXDxCrb5Y+BPgReSPD+M/SWwD3gsyd3A68Dtw7qngJ3AUeCXwGfXs2BJ0vKWDfeq+jcgS6y+aZHtC7hnzLokSWPwG6qS1JDhLkkNGe6S1JDhLkkNGe6S1JDhLkkNGe6S1JDhLkkNGe6S1NBKLj8g/b/ZvU9OZL/H9t06kf1K08ozd0lqyHCXpIYMd0lqyHCXpIYMd0lqyHCXpIYMd0lqyHCXpIYMd0lqyHCXpIYMd0lqyHCXpIYMd0lqyHCXpIYMd0lqyHCXpIYMd0lqyHCXpIYMd0lqyHCXpIYMd0lqyHCXpIYMd0lqyHCXpIYMd0lqyHCXpIaWDfckDyU5neTFkbH7kpxI8vxw2zmy7gtJjiZ5NcmnN6pwSdLSVnLm/k3g5kXG/7qqrh1uTwEk2Q7cAfzB8Jy/TXLRehUrSVqZZcO9qp4B3lrhz9sFPFpVv6qqnwFHgevGqE+StAbjzLnfm+TIMG1z+TB2JfDGyDbHh7H3SbInyeEkh8+cOTNGGZKkhdYa7g8AHweuBU4CX13tD6iq/VU1V1VzMzMzayxDkrSYNYV7VZ2qqveq6tfAN/jN1MsJYOvIplcNY5Kk82hN4Z5ky8jDzwBnP0lzELgjySVJrgG2AT8Yr0RJ0mpdvNwGSR4BbgQ2JTkOfAm4Mcm1QAHHgD8DqKqXkjwGvAy8C9xTVe9tSOWSpCUtG+5Vdeciww+eY/svA18epyhJ0nj8hqokNWS4S1JDhrskNWS4S1JDhrskNWS4S1JDhrskNWS4S1JDhrskNWS4S1JDhrskNWS4S1JDhrskNWS4S1JDhrskNWS4S1JDhrskNWS4S1JDhrskNWS4S1JDhrskNWS4S1JDF0+6gHHN7n1y0iVI0gXHM3dJashwl6SGDHdJashwl6SGDHdJashwl6SGDHdJashwl6SGDHdJashwl6SGpv7yA/pgmORlJo7tu3Vi+5bWyjN3SWrIcJekhpYN9yQPJTmd5MWRsSuSPJ3kteH+8mE8Se5PcjTJkSQ7NrJ4SdLiVnLm/k3g5gVje4FDVbUNODQ8BrgF2Dbc9gAPrE+ZkqTVWDbcq+oZ4K0Fw7uAA8PyAeC2kfGHa96zwGVJtqxTrZKkFVrrnPvmqjo5LL8JbB6WrwTeGNnu+DD2Pkn2JDmc5PCZM2fWWIYkaTFjv6FaVQXUGp63v6rmqmpuZmZm3DIkSSPWGu6nzk63DPenh/ETwNaR7a4axiRJ59Faw/0gsHtY3g08MTJ+1/CpmeuBt0embyRJ58my31BN8ghwI7ApyXHgS8A+4LEkdwOvA7cPmz8F7ASOAr8EPrsBNUuSlrFsuFfVnUusummRbQu4Z9yiJEnj8RuqktSQ4S5JDRnuktSQ4S5JDRnuktSQ4S5JDRnuktSQ4S5JDRnuktSQ4S5JDRnuktSQ4S5JDRnuktSQ4S5JDRnuktSQ4S5JDRnuktSQ4S5JDRnuktSQ4S5JDRnuktSQ4S5JDRnuktSQ4S5JDRnuktSQ4S5JDRnuktSQ4S5JDV086QKkC93s3icnst9j+26dyH7Vg2fuktSQ4S5JDRnuktSQ4S5JDRnuktSQ4S5JDRnuktTQWJ9zT3IMeAd4D3i3quaSXAH8IzALHANur6r/Gq9MSdJqrMeZ+59U1bVVNTc83gscqqptwKHhsSTpPNqIaZldwIFh+QBw2wbsQ5J0DuOGewHfTfJckj3D2OaqOjksvwlsHnMfkqRVGvfaMjdU1Ykkvwc8neTfR1dWVSWpxZ44/DHYA3D11VePWYYkadRYZ+5VdWK4Pw08DlwHnEqyBWC4P73Ec/dX1VxVzc3MzIxThiRpgTWHe5KPJLn07DLwKeBF4CCwe9hsN/DEuEVKklZnnGmZzcDjSc7+nH+oqn9O8kPgsSR3A68Dt49fpiRpNdYc7lX1U+APFxn/T+CmcYqSJI3Hb6hKUkOGuyQ1ZLhLUkP+D1XpAuX/btU4PHOXpIYMd0lqyHCXpIYMd0lqyHCXpIYMd0lqyHCXpIYMd0lqyHCXpIYMd0lqyHCXpIYMd0lqyHCXpIYMd0lqyHCXpIYMd0lqyHCXpIYMd0lqyH+zJ+m3TOrf+4H/4m89eeYuSQ0Z7pLUkOEuSQ0Z7pLUkOEuSQ0Z7pLUkOEuSQ0Z7pLUkOEuSQ0Z7pLUkJcfkPSB1/GSC565S1JDhrskNWS4S1JDGzbnnuRm4G+Ai4C/q6p9G7UvST1Mcu67mw05c09yEfB14BZgO3Bnku0bsS9J0vtt1LTMdcDRqvppVf0P8Ciwa4P2JUlaYKOmZa4E3hh5fBz4o9ENkuwB9gwP/zvJq2vc1ybg52t87oXO3qZT19669gUT7C1fGevpv7/Uiol9zr2q9gP7x/05SQ5X1dw6lHTBsbfp1LW3rn1Bz942alrmBLB15PFVw5gk6TzYqHD/IbAtyTVJPgzcARzcoH1JkhbYkGmZqno3yb3AvzD/UciHquqljdgX6zC1cwGzt+nUtbeufUHD3lJVk65BkrTO/IaqJDVkuEtSQ1Md7kluTvJqkqNJ9k66nnElOZbkhSTPJzk8jF2R5Okkrw33l0+6zpVI8lCS00leHBlbtJfMu384jkeS7Jhc5ee2RF/3JTkxHLfnk+wcWfeFoa9Xk3x6MlWvTJKtSb6X5OUkLyX53DDe4bgt1VuLY7eoqprKG/Nv1P4H8DHgw8CPge2TrmvMno4BmxaM/RWwd1jeC3xl0nWusJdPAjuAF5frBdgJ/BMQ4Hrg+5Ouf5V93Qf8xSLbbh9+Ly8Brhl+Xy+adA/n6G0LsGNYvhT4ydBDh+O2VG8tjt1it2k+c/+gXOJgF3BgWD4A3Da5Ulauqp4B3lowvFQvu4CHa96zwGVJtpyXQldpib6Wsgt4tKp+VVU/A44y/3t7Qaqqk1X1o2H5HeAV5r9t3uG4LdXbUqbq2C1mmsN9sUscnOtgTYMCvpvkueHyDACbq+rksPwmsHkypa2LpXrpcCzvHaYmHhqZOpvavpLMAp8Avk+z47agN2h27M6a5nDv6Iaq2sH81TTvSfLJ0ZU1/3qxxWdXO/UCPAB8HLgWOAl8daLVjCnJR4FvA5+vql+Mrpv247ZIb62O3ahpDvd2lzioqhPD/WngceZfBp46+1J3uD89uQrHtlQvU30sq+pUVb1XVb8GvsFvXr5PXV9JPsR8+H2rqr4zDLc4bov11unYLTTN4d7qEgdJPpLk0rPLwKeAF5nvafew2W7giclUuC6W6uUgcNfw6YvrgbdHpgEueAvmmT/D/HGD+b7uSHJJkmuAbcAPznd9K5UkwIPAK1X1tZFVU3/cluqty7Fb1KTf0R3nxvy79T9h/p3sL066njF7+Rjz787/GHjpbD/A7wKHgNeAfwWumHStK+znEeZf5v4v8/OVdy/VC/Oftvj6cBxfAOYmXf8q+/r7oe4jzIfClpHtvzj09Spwy6TrX6a3G5ifcjkCPD/cdjY5bkv11uLYLXbz8gOS1NA0T8tIkpZguEtSQ4a7JDVkuEtSQ4a7JDVkuEtSQ4a7JDX0f2xvmIv5cIVrAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# plot the emissions\n", "import matplotlib.pyplot as plt\n", "\n", "plt.hist(my_tier.co2_emissions.value)" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## Analytical error propagation\n", "\n", "Another method that is implemented in the ipcc package is analytical error propgation. The keyword `analytical` can be used in this case. The implementation is done using the `uncertainties` Python package. Thus all results are of type `ufloat`." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "2023-12-19 15:35:43,008 - INFO - Incineration sequence started --->\n", "/Users/TN76JP/Documents/coderefinery/ipcc/ipcc/src/ipcc/_sequence.py:87: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`\n", " new_c = conc_df.loc[c][j]\n", "2023-12-19 15:35:43,010 - INFO - 'Coordinates (2010, 'DE')' has been replaced by '[2006, 'Western Europe']' during reading parameter table 'msw_gen_rate'\n", "2023-12-19 15:35:43,011 - INFO - 'Coordinates (2010, 'DE')' has been replaced by '[2006, 'Western Europe']' during reading parameter table 'msw_frac_to_incin'\n", "2023-12-19 15:35:43,013 - INFO - 'Coordinates (2010, 'DE', 'msw_paper')' has been replaced by '[2006, 'DE', 'msw_paper']' during reading parameter table 'msw_type_frac'\n", "2023-12-19 15:35:43,015 - INFO - 'Coordinates (2010, 'DE', 'inc_unspecified')' has been replaced by '[2006, 'World', 'inc_unspecified']' during reading parameter table 'incintype_frac'\n", "2023-12-19 15:35:43,017 - INFO - 'Coordinates (2010, 'DE', 'msw_paper')' has been replaced by '[2006, 'World', 'msw_paper']' during reading parameter table 'dm'\n", "2023-12-19 15:35:43,019 - INFO - 'Coordinates (2010, 'DE', 'msw_paper')' has been replaced by '[2006, 'World', 'msw_paper']' during reading parameter table 'cf'\n", "2023-12-19 15:35:43,020 - INFO - 'Coordinates ('DE', 'msw_paper', 'inc_unspecified')' has been replaced by '['World', 'msw_paper', 'inc_unspecified']' during reading parameter table 'fcf'\n", "2023-12-19 15:35:43,020 - INFO - 'Coordinates ('DE', 'msw_paper', 'inc_unspecified')' has been replaced by '['World', 'msw_paper', 'inc_unspecified']' during reading parameter table 'of'\n", "2023-12-19 15:35:43,021 - INFO - ---> Incineration sequence finalized.\n" ] }, { "data": { "text/plain": [ "Step(position=12, year=2010, unit='Gg/year', value=88.4667825477414+/-52.63985799388056, type='elementary')" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "final_step = my_ipcc.waste.incineration.sequence.tier1_co2(year=2010, region=\"DE\", wastetype= \"msw_paper\", incintype= \"inc_unspecified\", uncertainty=\"analytical\").co2_emissions\n", "final_step" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Plot the results" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[]" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD4CAYAAADlwTGnAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAAsTAAALEwEAmpwYAAAwO0lEQVR4nO3deXxV9Z3/8dcnO0tIIAlbFhIggAFkC2EVF1QQlbigggvoYHUs1jqO7eB0potT29pOte24VUQFtSBubawKriDIEoIsEiCQlRAICQlZSMj+/f1xj/7SNJAL3OTc5fN8PPLw5pzvufd9vOF+7jnne75fMcaglFLK9/jZHUAppZQ9tAAopZSP0gKglFI+SguAUkr5KC0ASinlowLsDnAuIiMjTXx8vN0xlFLKY+zYseOEMSaqvXUeVQDi4+PJyMiwO4ZSSnkMESk40zo9BaSUUj5KC4BSSvkoLQBKKeWjtAAopZSP0gKglFI+SguAUkr5KC0ASinlozzqPgCl3MXJmgYyCk6y72gVzS0tAAT4+zEquhcT4voQ1j3Q5oRKdUwLgFJOqqlvYlX6YdZkFHLw+Knvlos4/vvt1BoiMLxfKAtS4rhtYiwhgf42pFWqY+JJE8IkJycbvRNYdbWqukZe3pTHq5vzqahtZGJ8by4b3peJ8X24OCbsuw/4usZmdh6uYHt+OZ8fKGFXYQWRPYO4Z1oC90yLp3uQft9SXU9EdhhjkttdpwVAqTPbllvGI2t2U1RxmquS+vH9y4YwLq63U9um55Xz7BfZbDhYSkJkD56+bSxjY8M7N7BSbWgBUOocNTS18IdPD/L8hhzi+nTnqVvHMmGQcx/8bW3OOcGja3ZzvLqeh2cm8sBlQwjw1/4XqmucrQDoX6FSbZyqb+LuV9J5bn0Ot06I5YOHLjnvD3+AqUMi+ejhGVw7egC//+Qg97+2g7rGZhcmVur8aAFQqpWTNQ3csWwr2/LK+d9bxvDkvIvpGXzh5+7DugXypwXj+J/UkXyeVcLCl9Opqmt0QWKlzp8WAKUsxZV13PLnLewvrubPd05g3oQYl7/GXVPi+eP8cXxdcJIFL27lxKl6l7+GUs7SAqAUUFnbyJ3Lt1FcWceKe1K4Mqlfp73W3DEDWbYwmZzSUyx6OZ2a+qZOey2lzkYLgPJ59U3NfO+1DA6X1fLSomSmDIno9Ne8fERfnr9jAvuPVbHkL1/T1NzS6a+pVFtaAJRPa2kxPPrWHtLzyvndLRczeXDnf/h/6/IRffnlDaNZn1XKf/11L57UI095B70zRfm0pz45yPu7j/Ifs0eQOja6y1//9klxFFXU8uwXOcRH9uBfLx3S5RmU79IjAOWzvsgq4ZkvsrktOZZ/vXSwbTkevXo41148gN+uPUB6XrltOZTvcaoAiMhsEckSkWwRWdrO+mARedNav01E4lute8xaniUis6xlw0VkV6ufKhF52FU7pVRHiivr+Pc1uxnRP5RfpI5Evh3QxwYiwm9uGk1cn+48tGon5TUNtmVRvqXDAiAi/sCzwDVAErBARJLaNFsMnDTGDAWeBp60tk0C5gMjgdnAcyLib4zJMsaMNcaMBSYAtcB7rtklpc6uqbmFh1btpK6xmWfvGO8Wg7WFhgTyzO3jKa9p4JE1u2hp0esBqvM5cwSQAmQbY3KNMQ3AaiC1TZtUYIX1+G1gpji+UqUCq40x9caYPCDber7WZgI5xpiC890Jpc7Fnz47RHp+OU/cOIohUT3tjvOdUdFh/Pd1F7E+q5RlG3PtjqN8gDMFIBoobPX7EWtZu22MMU1AJRDh5LbzgVVnenERuU9EMkQko7S01Im4Sp3ZN0cqeXZ9DjePj+HGca6/0etC3Tl5ELNH9uf3Hx8ku6Ta7jjKy9l6EVhEgoC5wFtnamOMedEYk2yMSY6Kiuq6cMrrNDS18KO3dxPRI4ifXt/2LKZ7EBH+54ZRdA/258dv76FZTwWpTuRMASgCYlv9HmMta7eNiAQAYUCZE9teA3xtjDl+brGVOnfPr8/hQHE1T9w4mrBu7jtjV1RoMD+7PomvD1fw6uZ8u+MoL+ZMAdgOJIpIgvWNfT6Q1qZNGrDIejwP+Nw47mpJA+ZbvYQSgEQgvdV2CzjL6R+lXCWruJpnvjjE3DEDuaoTh3lwlRvGRnPFiL78bt0BCspq7I6jvFSHBcA6p/8gsA7YD6wxxmSKyOMiMtdqthyIEJFs4BFgqbVtJrAG2AesBZYYY5oBRKQHcBXwrmt3Sal/1NJi+I939hAaEsjP3PTUT1siwhM3jiLQz4//fO8bvUtYdQqn7gQ2xnwIfNhm2U9bPa4DbjnDtk8AT7SzvAbHhWKlOtXbXx9hV2EFv79lDBE9g+2O47QBYd340ezh/PRvmXy0t5g5owfYHUl5Gb0TWHm1qrpGfrv2AOPjwrlpfNcP9XChbk+JY0T/UJ74YD+nG3QSGeVaWgCUV/vDJ4coq2ng8dRRtt7te74C/P34+dyRFFWc5vkNOXbHUV5GC4DyWgePV7NiSz7zJ8YxKjrM7jjnbfLgCK4fM5AXNuRQWF5rdxzlRbQAKK9kjOEX72fSMziAH80abnecC/afc0bgL8IvP9hndxTlRbQAKK+04WApX2WX8fCVifTpEWR3nAs2IKwb379sCOsyj5ORryOGKtfQAqC8TnOL4TcfHSCuT3fumDTI7jgus/iSBPqGBvPrjw5ot1DlEloAlNf5264iDhRX8+is4QQFeM+fePegAB6+chg7Ck7yyT69eV5dOO/516EUUNfYzO8/Psio6F5c54X95m9NjmFwZA9+uy5L5xFWF0wLgPIqr28toKjiNEtnX4Sfn+d1++xIgL8fP549nOySU7zz9RG74ygPpwVAeY3qukae/SKbSxIjmZ4YaXecTjNrZH/GxYXz9CeHqGvUm8PU+dMCoLzGis35nKxt5NGrPb/b59mICD+6ejjFVXWsTj9sdxzlwbQAKK9QXdfIso15zBzRlzGx4XbH6XRThkSQktCH59bn6FGAOm9aAJRXePWrfCpPN/LwlcPsjtIlRIR/u3IYJdX1rNKjAHWetAAoj1dV18iyjblceVFfRsd47pAP52rKkAgm6VGAugBaAJTHe2VTPlV1TT7z7b+1f7tqGKXV9byxTY8C1LnTAqA8WlVdI8s35XJVUj+PHvDtfE0eHMGUwRE8r0cB6jxoAVAe7fWtBVTVNfHQFYl2R7HNQzMTOXGqnrcyCu2OojyMFgDlseoam3l5Ux4zhkX51Ln/tiYP7sP4uHBe2JBLo94drM6BUwVARGaLSJaIZIvI0nbWB4vIm9b6bSIS32rdY9byLBGZ1Wp5uIi8LSIHRGS/iExxyR4pn7Emo5ATpxpYctkQu6PYSkRYcvlQiipO8/7uo3bHUR6kwwIgIv7As8A1QBKwQETazqy9GDhpjBkKPA08aW2bBMwHRgKzgees5wP4I7DWGDMCGINjwnmlnNLY3MKfN+QyYVBvUhL62B3HdleM6MuI/qE8tz6HlhYdKVQ5x5kjgBQg2xiTa4xpAFYDqW3apAIrrMdvAzPFMf9eKrDaGFNvjMkDsoEUEQkDZgDLAYwxDcaYigveG+Uz0nYdpajiNEsuH+KRUz26mojwwGVDyC45xcc6UqhykjMFIBpofXXpiLWs3TbGmCagEog4y7YJQCnwiojsFJGXRKRHey8uIveJSIaIZJSWljoRV3m7lhbD8xtyGNE/lMuH97U7jtu4dvQA4vp057n12TpfgHKKXReBA4DxwPPGmHFADfBP1xYAjDEvGmOSjTHJUVFRXZlRualP9x8nu+QUD1ym3/5bC/D34/5LB7PnSCVbcsvsjqM8gDMFoAiIbfV7jLWs3TYiEgCEAWVn2fYIcMQYs81a/jaOgqBUh5ZtzCU6vBvXeuF4/xfq5vExRPYM4sUvc+2OojyAMwVgO5AoIgkiEoTjom5amzZpwCLr8Tzgc+M4Bk0D5lu9hBKARCDdGFMMFIrIt8M2zgR0tmvVoa8Pn2R7/kkWT08gwF97MbcVEujPwinxrM8qJau42u44ys11+C/IOqf/ILAOR0+dNcaYTBF5XETmWs2WAxEikg08gnU6xxiTCazB8eG+FlhijPn2dsUfAG+IyB5gLPArl+2V8lrLvsylV0gAt02M7bixj7pr8iBCAv14aaMeBaizC3CmkTHmQ+DDNst+2upxHXDLGbZ9AniineW7gORzyKp8XEFZDWszi3ng0iH0CHbqT9cn9e4RxK3JsaxKP8yjs4bTr1eI3ZGUm9JjaOUxXtqYR6CfH3dPjbc7ittbPD2B5hbDq5vz7Y6i3JgWAOURymsaeGtHITeMG0hf/UbboUERPZg9qj+vby3gVH2T3XGUm9ICoDzCX7YVUNfYwr2XDLY7ise495LBVNc18bYOEqfOQAuAcnsNTS2s3FLAJYmRDOsXanccjzE+rjfj4sJ5ZXO+Dg+h2qUFQLm9D745Skl1Pf8yPcHuKB5n8fQECspq+exAid1RlBvSAqDcmjGG5ZvyGBLVg0sT9U7wczV7ZH8GhoXw8qY8u6MoN6QFQLm17fkn2VtUxT3TEvDz02EfzlWAvx+LpsazJbeMzKOVdsdRbkYLgHJrL2/KI6xbIDePj7E7iseaPzGOboH+vPJVvt1RlJvRAqDcVmF5LR/vK+b2SXF0C/LveAPVrrDugdySHEParqOUVtfbHUe5ES0Aym2t3JKPiLBwyiC7o3i8u6fG09Dcwl+2HbY7inIjWgCUW6ptaOLN7YXMHtWfAWHd7I7j8QZH9eSy4VG8sa2AhiadN1g5aAFQbum9nUVU1TVxjw774DJ3T42npLqej/YeszuKchNaAJTbMcawYnM+o6J7MWFQb7vjeI0ZiVEkRPbQ8YHUd7QAKLezJaeMg8dPcffUBJ3xy4X8/IRFUwax83AFuwsr7I6j3IAWAOV2XtmcT58eQVx3sc745Wo3T4ihZ3AAK/QoQKEFQLmZwvJaPt1/nNtT4ggJ1K6frhYaEsi8CTG8v0e7hCotAMrNvL61AD8R7pgcZ3cUr7VwyiAamw2r07VLqK9zqgCIyGwRyRKRbBFZ2s76YBF501q/TUTiW617zFqeJSKzWi3PF5FvRGSXiGS4ZG+URzvd0Mzq7YXMHqldPzvT4KiezBgWxevbCmhs1i6hvqzDAiAi/sCzwDVAErBARJLaNFsMnDTGDAWeBp60tk3CMYn8SGA28Jz1fN+63Bgz1hijU0Mq0nYXUXm6UW/86gKLpgzieFU9H2cetzuKspEzRwApQLYxJtcY0wCsBlLbtEkFVliP3wZmiqP7Riqw2hhTb4zJA7Kt51PqHzi6fhYwon8oKQl97I7j9S4b3pfYPt1YsSXf7ijKRs4UgGig9ZRCR6xl7bYxxjQBlUBEB9sa4GMR2SEi953pxUXkPhHJEJGM0tJSJ+IqT7Sj4CT7jlWxcEq8dv3sAv5+wl2TB5GeV87+Y1V2x1E2sfMi8HRjzHgcp5aWiMiM9hoZY140xiQbY5KjonQ8eG+1YksBoSEB3DBuoN1RfMatybEEB/ixckuB3VGUTZwpAEVAbKvfY6xl7bYRkQAgDCg727bGmG//WwK8h54a8lklVXV89M0xbk2OpXtQgN1xfEZ49yBuGBvNX3cWUVnbaHccZQNnCsB2IFFEEkQkCMdF3bQ2bdKARdbjecDnxhhjLZ9v9RJKABKBdBHpISKhACLSA7ga2Hvhu6M80V/SD9PUYrhrsl787WoLpw7idGMzb+3QieN9UYcFwDqn/yCwDtgPrDHGZIrI4yIy12q2HIgQkWzgEWCptW0msAbYB6wFlhhjmoF+wCYR2Q2kAx8YY9a6dteUJ2i0hii+dFgU8ZE97I7jc0YODGPCoN68vrVAJ473QU4dbxtjPgQ+bLPsp60e1wG3nGHbJ4An2izLBcaca1jlfdZlFlNSXc+vb9Jv/3ZZOGUQP1y9iy8PlXLZ8L52x1FdSO8EVrZauaWAmN7d9IPHRteMGkBkz2Be04vBPkcLgLLNgeIq0vPKuXPyIPx1wnfbBAX4sSAlls+zSigsr7U7jupCWgCUbV7bUkBwgB+3Jcd23Fh1qtsnxeEnwuvb9CjAl2gBULaoqmvkvZ1FXD9mIL17BNkdx+cNCOvGVRf1Y832Quoam+2Oo7qIFgBli3d3HKG2oVnH/XEjC6cO4mRtI+/vPmp3FNVFtACoLmeMYeXWAsbEhnNxTLjdcZRlyuAIEvv25LWtehrIV2gBUF3uq+wycktrWKg3frkVEeGuKYPYc6SSXTplpE/QAqC63Motjikfr9UpH93OjeOi6RHkz0odJdQnaAFQXaqo4jSf7j/ObRNjdcpHNxQaEshN42P4+55jlJ3SKSO9nRYA1aXesM4v3zFJp3x0VwunDKKhqYU3M3R8IG+nBUB1mfqmZt7cXsgVI/oR07u73XHUGST2C2XK4Aje2HqYZh0fyKtpAVBd5sNvjlFW08CiqXrx190tnDKIoorTfH6gxO4oqhNpAVBdZuWWAgZH9mDakEi7o6gOXJXUj/69QvRisJfTAqC6xJ4jFew8XMFdUwbhp+P+uL0Afz/umBTHxkMnyCk9ZXcc1Um0AKgusXJLAd2D/Ll5QozdUZST5qfEEegvOkqoF9MCoDpdeU0DabuPctP4aHqFBNodRzkpKjSYa0cP4J0dRzhV32R3HNUJtACoTvfm9kIamlpYOCXe7ijqHC2cGk91fRPv7Ww7DbjyBk4VABGZLSJZIpItIkvbWR8sIm9a67eJSHyrdY9Zy7NEZFab7fxFZKeI/P2C90S5peYWw+tbC5gyOIJh/ULtjqPO0bjYcEZHh7Fycz6Oab6VN+mwAIiIP/AscA2QBCwQkaQ2zRYDJ40xQ4GngSetbZNwTCI/EpgNPGc937d+iGOeYeWlPj9QQlHFae366aFEhIVTBnGo5BRbcsvsjqNczJkjgBQg2xiTa4xpAFYDqW3apAIrrMdvAzNFRKzlq40x9caYPCDbej5EJAa4FnjpwndDuasVm/MZEBbClRf1szuKOk/XjxlI7+6BrNicb3cU5WLOFIBooPU94UesZe22McY0AZVARAfb/gH4MdByrqGVZzh0vJpN2Se4c/IgAvz1cpOnCgn0Z35KHJ/sO86RkzplpDex5V+liFwHlBhjdjjR9j4RyRCRjNLS0i5Ip1xlxZZ8a75ZHffH091pDd2tcwV4F2cKQBHQetLWGGtZu21EJAAIA8rOsu00YK6I5OM4pXSFiLze3osbY140xiQbY5KjoqKciKvcQeXpRt79uoi5YwbSR6d89HjR4d2YNbI/b+qUkV7FmQKwHUgUkQQRCcJxUTetTZs0YJH1eB7wuXF0GUgD5lu9hBKARCDdGPOYMSbGGBNvPd/nxpg7XbA/yk28lVFIbUMzd0+NtzuKcpFFU+OpqG3kb7u0S6i36LAAWOf0HwTW4eixs8YYkykij4vIXKvZciBCRLKBR4Cl1raZwBpgH7AWWGKM0a8PXq65xbBySwHJg3ozKjrM7jjKRSYl9GFE/1Be+Uq7hHqLAGcaGWM+BD5ss+ynrR7XAbecYdsngCfO8tzrgfXO5FCeYX1WCYfLa/nRrOF2R1EuJCLcPTWepe9+Q3peOZMGR9gdSV0g7ZqhXO7Vzfn06xXM7FH97Y6iXCx1bDTh3QN5VbuEegUtAMqlDh6vZuOhEyycEk+gdv30Ot2C/Jk/MY51mcXaJdQL6L9Q5VKvfJVPsHb99GoLpwxCRFipo4R6PC0AymVO1jTw3s4j3DguWrt+erGB4d24ZlR/VqUfpkZHCfVoWgCUy6zafpi6xhbumZZgdxTVye6ZlkB1XRPvfH3E7ijqAmgBUC7R2NzCys0FTB8ayfD+OuqntxsfF86Y2HBe/SqfFp043mNpAVAusXZvMcVVddwzLd7uKKoLiAj/Mi2e3BM1bDioQ7R4Ki0A6oIZY1i+KY/4iO5cPryv3XFUF5kzegD9egWzfFOe3VHUedICoC7YjoKT7Cqs4F+mJ+iE7z4k0N+PhVPi2ZR9gv3HquyOo86DFgB1wV7amEdYt0Dm6YTvPueOSXF0C/TnpY16FOCJtACoC1JQVsO6fcXcOTmO7kFOjSyivEh49yBuTY4hbXcRx6vq7I6jzpEWAHVBXt6UR4Cf6ITvPuyeaQk0tRhWbsm3O4o6R1oA1HmrqG1gTcYR5o6Jpl+vELvjKJvER/bg6qR+vL71MLUNemOYJ9ECoM7bG9sOc7qxmXsv0Ru/fN33LhlM5elG3t6hN4Z5Ei0A6rzUNzWzYnM+04dGctGAXnbHUTabMKg3Y2LDWb4pj2a9McxjaAFQ5+WvO4soqa7nvhmD7Y6i3ICIcP+MwRSU1bJ2b7HdcZSTtACoc9bSYvjzl7kkDejFJYmRdsdRbmLWyP7ER3TnhQ05OmOYh9ACoM7Zp/uPk1taw/2XDkZEb/xSDv5+wvdmDOabokq25JTZHUc5wakCICKzRSRLRLJFZGk764NF5E1r/TYRiW+17jFreZaIzLKWhYhIuojsFpFMEfmFy/ZIdbo/f5lLTO9uXDt6gN1RlJu5eXwMkT2DeOHLXLujKCd0WABExB94FrgGSAIWiEhSm2aLgZPGmKHA08CT1rZJwHxgJDAbeM56vnrgCmPMGGAsMFtEJrtkj1Sn2p5fzo6Ck3zvksEE6Ixfqo2QQH/umZbAlwdLyTxaaXcc1QFn/gWnANnGmFxjTAOwGkht0yYVWGE9fhuYKY5zA6nAamNMvTEmD8gGUozDKat9oPWjJw09wJ835NCnRxC3JsfaHUW5qTsnDaJHkD9/3qBHAe7OmQIQDRS2+v2ItazdNsaYJqASiDjbtiLiLyK7gBLgE2PMtvZeXETuE5EMEckoLdVhZ+10oLiKT/eXsGhKPN2C/O2Oo9xUWPdAbp8UxwffHONwmc4b7M5sO4Y3xjQbY8YCMUCKiIw6Q7sXjTHJxpjkqKioLs2o/tGzX+TQI8ifRVMH2R1Fubl7LxmMvwjPb8ixO4o6C2cKQBHQ+ng/xlrWbhsRCQDCgDJntjXGVABf4LhGoNxU3okaPthzlLumxBPeXef7VWfXr1cIt06M4Z0dRyiu1EHi3JUzBWA7kCgiCSIShOOiblqbNmnAIuvxPOBz4+gInAbMt3oJJQCJQLqIRIlIOICIdAOuAg5c8N6oTvP8+mwC/f1YPF2HfVDOuX/GEJqN4UXtEeS2OiwA1jn9B4F1wH5gjTEmU0QeF5G5VrPlQISIZAOPAEutbTOBNcA+YC2wxBjTDAwAvhCRPTgKzCfGmL+7dteUqxRVnObdr4tYkBJHVGiw3XGUh4jt050bxkbzl/QCyk7V2x1HtUM86Y695ORkk5GRYXcMn/Ozv+3lL+mH2fCjyxkY3s3uOMqD5JSe4sqnNvD9y4bwo1kj7I7jk0RkhzEmub112pFbnVVJVR2rtxdy07gY/fBX52xIVE/mjBrAys0FVNQ22B1HtaEFQJ3V8xtyaGoxLLl8qN1RlIf6wcyhVNc38bJOHu92tACoMzpeVccb2w5z8/ho4iK62x1HeagR/XsxZ3R/Xv4qX48C3IwWAHVGz6/PoaXF8ODliXZHUR7uhzOHcaq+SSePdzNaAFS7iivr+Ev6YW4eH6Pf/tUFG94/lGtHD+DVzfmcrNGjAHehBUC16/n12Y5v/1fouX/lGg/NTKSmoYmXNul9Ae5CC4D6J0crTrMqvZB5E2KI7aPf/pVrDO8fypzRA3j1q3y9L8BNaAFQ/+SPnx4C0G//yuX+7cpETjc289x6HSPIHWgBUP8gp/QUb+0o5PZJccT01m//yrWG9g3lpvExvLa1gKMVp+2O4/O0AKh/8NTHBwkJ9Ndv/6rTPHxlIpj/f6Sp7KMFQH1nb1ElH3xzjMXTE4jsqWP+qM4R07s7t0+K460dheSUnup4A9VptACo7/x2XRbh3QP53ozBdkdRXu7BK4YSEujPUx8ftDuKT9MCoADYnHOCLw+W8v3LhtArJNDuOMrLRfYM5t7pCXzwzTF2F1bYHcdnaQFQtLQYnvhgP9Hh3Vg4Jd7uOMpHfG/GYCJ7BvHEB/vxpFGJvYkWAMW7O4vIPFrFj2cPJyRQ5/pVXSM0JJCHrxxGen456zKP2x3HJ2kB8HGnG5r533VZXBwTxvUXD7Q7jvIx8yfGMrRvT37z0X4amlrsjuNztAD4uJc25lJcVcdP5lyEn5/YHUf5mAB/P/5zzgjyy2p5Y1uB3XF8jlMFQERmi0iWiGSLyNJ21geLyJvW+m0iEt9q3WPW8iwRmWUtixWRL0Rkn4hkisgPXbZHymkl1XW8sCGHq5P6MWlwhN1xlI+6fHhfpg2N4I+fHaKyttHuOD6lwwIgIv7As8A1QBKwQESS2jRbDJw0xgwFngaetLZNwjGJ/EhgNvCc9XxNwL8bY5KAycCSdp5TdbInP8qiobmFx+ZcZHcU5cNEhJ/MSaLqdCNPf6rdQruSM0cAKUC2MSbXGNMArAZS27RJBVZYj98GZoqIWMtXG2PqjTF5QDaQYow5Zoz5GsAYU41jsvnoC98d5awdBSd55+sjLJ4+mITIHnbHUT4uaWAv7pg0iJVb8tl/rMruOD7DmQIQDRS2+v0I//xh/V0bY0wTUAlEOLOtdbpoHLCtvRcXkftEJENEMkpLS52IqzrS3GL4Wdpe+vUK5gc65INyE/9+9TDCugXys7RM7RbaRWy9CCwiPYF3gIeNMe2WfWPMi8aYZGNMclRUVNcG9FJvbi9kb1EV/znnInoEB9gdRykAwrsH8eis4aTnlfP+nmN2x/EJzhSAIiC21e8x1rJ224hIABAGlJ1tWxEJxPHh/4Yx5t3zCa/OXUVtA79bd4CUhD7MHaPdPpV7mT8xjlHRvfjVB/upqW+yO47Xc6YAbAcSRSRBRIJwXNRNa9MmDVhkPZ4HfG4cx3BpwHyrl1ACkAikW9cHlgP7jTFPuWJHlHN+89EBquqa+MXckTjeBqXch7+f8Iu5oyiuquPpT/SCcGfrsABY5/QfBNbhuFi7xhiTKSKPi8hcq9lyIEJEsoFHgKXWtpnAGmAfsBZYYoxpBqYBdwFXiMgu62eOi/dNtbEtt4zV2wu5d3oCFw3oZXccpdo1YVBvbp8Ux8tf5fHNkUq743g18aSLLcnJySYjI8PuGB6pvqmZa/64kYamFj7+txl0D9Jz/8p9VZ5u5MqnNtCvVzB//f40Avz1ntXzJSI7jDHJ7a3T/6s+4vn1OeSW1vDLG0bph79ye2HdAvn59SPZW1TFq5vz7Y7jtbQA+IDskmqe+yKHuWMGctnwvnbHUcopc0b3Z+aIvvz+44MUltfaHccraQHwck3NLTz61h66B/vz39fpzdbKc4gIj98wCj+Bpe/uoaXFc05XewotAF7uxY257Cqs4PHUUUSF6jSPyrNEh3fjv65L4qvsMh0srhNoAfBiWcXV/OGTQ8wZ3Z/rLx5gdxylzsv8ibFcOiyKX314gIKyGrvjeBUtAF6qsbmFR9bsIjQkgP9JHaV9/pXHEhF+c/NoAvyFH72lp4JcSQuAl/rTZ4fIPFrFr24aTURPPfWjPNuAsG78/PqRpOeXs2xjrt1xvIYWAC+0JaeMZ77IZt6EGGaN7G93HKVc4qbx0Vwzqj+/W5elE8m7iBYAL1Ne08DDb+4kIbIHv5g70u44SrmMiPCbmy6mX68QfrBqJ9V1OnnMhdIC4EWMMfz47d2crGnkT/PH6UifyuuEdQ/kj/PHcuRkLf/11706bPQF0gLgRV75Kp9P95ew9JoRjIoOszuOUp0iOb4PD185jL/tOspbO47YHcejaQHwEul55fzqw/1ceVE/7pkWb3ccpTrVksuHMnVIBP/1173sLdIB486XFgAvcLyqju+/8TWxfbrz1G1jtMun8nr+fsL/LRhHZI8g7n9tB+U1DXZH8khaADxcQ1MLD7y+g9qGJv581wR6hQTaHUmpLhHRM5gX7ppA6al6Hlq1k2a9P+CcaQHwYMYYfpaWydeHK/jdvDEM6xdqdySlutTFMeH8MnUUm7JP8JuP9tsdx+NoNxEPtmxjLqvSD/PAZUO4Vod6UD7q1omx7D1aybKNeSRE9uT2SXF2R/IYWgA81Nq9xfz6owNcO3oAP7p6uN1xlLLVT69L4nB5Lf/9t73E9O7GjGFRdkfyCHoKyAPtLqzg4Td3MjY2nN/fOgY/P73oq3xbgL8fz9w+nsS+PVnyxtdkFVfbHckjOFUARGS2iGSJSLaILG1nfbCIvGmt3yYi8a3WPWYtzxKRWa2WvywiJSKy1yV74iOyS05xz6vbiewZzLKFyYQE+tsdSSm30DM4gJfvnki3IH8WvrxNJ5FxQocFQET8gWeBa4AkYIGItJ1ZZDFw0hgzFHgaeNLaNgmYD4wEZgPPWc8H8Kq1TDnpyMla7lq+DT8RXls8iUgd5E2pfzAwvBsrF6dQ19jCncu3UVJVZ3ckt+bMEUAKkG2MyTXGNACrgdQ2bVKBFdbjt4GZ4uiMngqsNsbUG2PygGzr+TDGfAmUu2AffEJpdT13vrSNmvomXlucQkJkD7sjKeWWRvTvxSv3TKS0up67lqdTUav3CJyJMwUgGihs9fsRa1m7bYwxTUAlEOHktmclIveJSIaIZJSWlp7Lpl7j2w//41X1vHJPChcN6GV3JKXc2vi43ixbmEzeiRotAmfh9heBjTEvGmOSjTHJUVG+d2W/uLKO217cwuHyWpYvSmbCoN52R1LKI0wbGskLd40nq7iaBcu2UXaq3u5IbseZAlAExLb6PcZa1m4bEQkAwoAyJ7dVZ1BUcZrbXtzC8co6VvxLClOHRtodSSmPcsWIfry0KJnc0lMsWLaVkmq9JtCaMwVgO5AoIgkiEoTjom5amzZpwCLr8Tzgc+MYpzUNmG/1EkoAEoF010T3bgePV3PL85spr2ngtXsnkZLQx+5ISnmkGcOieOWeiRSWn+bWF7bovMKtdFgArHP6DwLrgP3AGmNMpog8LiJzrWbLgQgRyQYeAZZa22YCa4B9wFpgiTGmGUBEVgFbgOEickREFrt21zzX1twy5j2/mcYWw6rvTWZ8nJ72UepCTB0Syev3TqLidCM3PbeZXTqjGADiSRMqJCcnm4yMDLtjdKq03Ud5dM1uYvt049V7Uojt093uSEp5jZzSU9z9Sjql1fX834LxXJXUz+5InU5Edhhjkttb5/YXgX1Fc4vht2sP8NCqnYyJDeOdB6bqh79SLjYkqifvPjCNYf1Cue+1DJ75/BAtPjyKqBYAN1BR28A9r27nufU5LEiJ4/V7JxHePcjuWEp5pajQYN68bwqpYwbyvx8f5IE3dnCqvsnuWLbQAmCznYdPcv0zm9iaU8avbxrNr28aTXCADu+gVGfqFuTP07eN5b+vS+LT/SXMfWYTmUd9b2YxLQA2aW4x/N9nh5j3whZaWmD1/ZNZkKLD2CrVVUSExdMTeH3xJGrqm7jx2c28tDHXp04JaQGwQd6JGha8uJXff3KQ6y4ewEcPX6I9fZSyyZQhEaz94QwuGx7FLz/Yz8KX0zly0jcGktNeQF2osbmFZRtz+cOnhwgO8OMXc0dy0/gYu2MppXDMsLcqvZAnPtiHAR69ejiLpsbj7+HDrZ+tF5AWgC6yLbeMn7+/j/3HqrhmVH9+MXckfXuF2B1LKdVGUcVpfvLeN6zPKmVMTBg/nzuScR58hK4FwEaHy2r59Uf7+WhvMQPCQvjZ9SOZPaq/3bGUUmdhjCFt91F++cF+SqvruWHsQH48ewQDw7vZHe2cna0A6JSQneR4VR3PfZHNqvRC/P2ER64axvcuGUy3IO3ho5S7ExFSx0Yz86J+PL8+m2Ub81ibWcxdkwdx/6VDvGYuDj0CcLFjladZ9mUer28roLnFMG98DA9flciAMM/75qCUcigsr+WpTw7yt11FBAf4s3DqIBZPS/CI07h6CqgL7C2qZNnGXD7YcwwD3DQumh9ckUhchN7Nq5S3yCk9xZ8+O0Ta7qME+vkxd+xA7r0kgRH93XeODi0AnaSmvom/7znKqvRCdhVW0DM4gNsmxnL31HgdxkEpL5Z3ooZXvsrjrYwjnG5sJiW+D7dNjGXO6AFud5pXC4ALNTa3sCn7BO/vPsq6vcXUNDQztG9P5k+M5daJsfQKCbQ1n1Kq61TUNrB6eyGr0w+TX1ZLaEgAc0YN4PoxA5k8uA8B/vbfaqUF4ALV1Dex8dAJPtt/nE/3H+dkbSOhIQHMHtmf2ybGMmFQbxxTICulfJExhq255azJKOTjTMcXw8iewVyV1JeZI/oxbWikbUcG2gvoHDU1t7D3aBVfZZ9gc84JtuedpKG5hdCQAK4Y0ZfrLh7IjGGROmaPUgpw9BqaMiSCKUMiqGts5osDJfx9zzHe332MVemFBAf4kZLQh2lDI5k2JJKkgb3c4gYznz8CMMZwvKqevUWV7D5SwY6Ck+wqrKC2oRmAEf1DmT40kpkX9SM5vjeBbnBIp5TyDA1NLWzLK+Oz/SV8lX2CQyWnAOgZHMC4uHDGx/VmbGw4I6N70Te0c3oU6SkgHB/0ZTUN5J+oIbvkFAePn+JQSTX7j1Vx4lQDAP5+wkUDQpkQ15sJ8X2YOiTCa/r7KqXsV1JVx+acMjIKytlRUEFWcRXfjj3XNzSYEQN6MaxvT4b1C2VI354kRPagd/fACzrFfMEFQERmA38E/IGXjDG/abM+GFgJTMAxGfxtxph8a91jwGKgGXjIGLPOmedsz/kUgOYWw43PfUVeaQ3Vrcb8Dgn0I7FvKMP7hzJqYC9GRYeRNLAX3YP0rJhSqmucqm9i39Eq9hZVsreokgPF1eSUnqK+qeW7Nr1CAhjRvxdv3j/5vArBBV0DEBF/4FngKuAIsF1E0owx+1o1WwycNMYMFZH5wJPAbSKShGMS+ZHAQOBTERlmbdPRc7qEv58wOLIH42LDiY/sQXxED4ZE9SSmdzf83OAcnFLKd/UMDiAloQ8pCX2+W9bcYigsryX3xCnyTtRSUFZDY3NLp3Q0cebrbgqQbYzJBRCR1UAqjonev5UK/Nx6/DbwjDjSpgKrjTH1QJ41aXyK1a6j53SZP8wf1xlPq5RSLufvJ44vq5E9Ov21nLmiGQ0Utvr9iLWs3TbGmCagEog4y7bOPCcAInKfiGSISEZpaakTcZVSSjnD7bu0GGNeNMYkG2OSo6Ki7I6jlFJew5kCUATEtvo9xlrWbhsRCQDCcFwMPtO2zjynUkqpTuRMAdgOJIpIgogE4biom9amTRqwyHo8D/jcOLoXpQHzRSRYRBKARCDdyedUSinViTq8CGyMaRKRB4F1OLpsvmyMyRSRx4EMY0wasBx4zbrIW47jAx2r3RocF3ebgCXGmGaA9p7T9bunlFLqTHzmRjCllPJFZ7sPwO0vAiullOocWgCUUspHedQpIBEpBQrszgFEAifsDnEBNL+9NL+9PD0/nNs+DDLGtNuH3qMKgLsQkYwznVPzBJrfXprfXp6eH1y3D3oKSCmlfJQWAKWU8lFaAM7Pi3YHuECa316a316enh9ctA96DUAppXyUHgEopZSP0gKglFI+SguAk0Tk5yJSJCK7rJ85rdY9JiLZIpIlIrPszHk2IjLbypgtIkvtzuMMEckXkW+s/+cZ1rI+IvKJiByy/tvb7pyticjLIlIiIntbLWs3szj8yXpP9ojIePuSf5e1vfwe8/cvIrEi8oWI7BORTBH5obXcI96Ds+R3/XtgjNEfJ35wzHj2aDvLk4DdQDCQAOQA/nbnbSenv5VtMBBkZU6yO5cTufOByDbLfgsstR4vBZ60O2ebfDOA8cDejjIDc4CPAAEmA9vcNL/H/P0DA4Dx1uNQ4KCV0yPeg7Pkd/l7oEcAF+67aS+NMXlA62kv3cl3U3saYxqAb6fh9ESpwArr8QrgBvui/DNjzJc4RsVt7UyZU4GVxmErEC4iA7ok6BmcIf+ZuN3fvzHmmDHma+txNbAfx4yDHvEenCX/mZz3e6AF4Nw8aB0ivtzqtIPT01vazFNytmWAj0Vkh4jcZy3rZ4w5Zj0uBvrZE+2cnCmzJ70vHvf3LyLxwDhgGx74HrTJDy5+D7QAtCIin4rI3nZ+UoHngSHAWOAY8Hs7s/qQ6caY8cA1wBIRmdF6pXEcA3tUX2ZPzIwH/v2LSE/gHeBhY0xV63We8B60k9/l70GHE8L4EmPMlc60E5FlwN+tXz1lektPyfkPjDFF1n9LROQ9HIe2x0VkgDHmmHWoXmJrSOecKbNHvC/GmOPfPvaEv38RCcTx4fmGMeZda7HHvAft5e+M90CPAJzU5pzgjcC3PSTONO2lu/G4aThFpIeIhH77GLgax//31lOQLgL+Zk/Cc3KmzGnAQqsnymSgstVpCrfhSX//IiI4Zincb4x5qtUqj3gPzpS/U94DO692e9IP8BrwDbDH+h8+oNW6n+C48p4FXGN31rPswxwcPQpygJ/YnceJvINx9G7YDWR+mxmIAD4DDgGfAn3sztom9yoch+iNOM7HLj5TZhw9T5613pNvgGQ3ze8xf//AdBynd/YAu6yfOZ7yHpwlv8vfAx0KQimlfJSeAlJKKR+lBUAppXyUFgCllPJRWgCUUspHaQFQSikfpQVAKaV8lBYApZTyUf8PTwib2UfLu+sAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "import scipy.stats as stats\n", "\n", "mu = final_step.value.n\n", "sigma = final_step.value.std_dev\n", "x = np.linspace(mu - 3*sigma, mu + 3*sigma, 100)\n", "plt.plot(x, stats.norm.pdf(x, mu, sigma))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## What can we learn?\n", "\n", "In this example, the values can be negative when using analytical error propagation. This is due to the high uncertainties in the data, using the tier 1 approach. However, it does not make sense to have negative CO2 emissions! Thus, Monte Carlo simulation should be used rather than analytical error propagation which is suitable for lower uncertainties." ] } ], "metadata": { "kernelspec": { "display_name": "ipcc", "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.10.4" }, "orig_nbformat": 4 }, "nbformat": 4, "nbformat_minor": 2 }