{ "cells": [ { "cell_type": "markdown", "id": "320a2808", "metadata": {}, "source": [ "# Chance-constrained PULPO under uncertainty: an open example\n", "\n", "This notebook runs PULPO's uncertainty method end to end on the bundled\n", "`soc_demo_database`: a toy ammonia system whose hydrogen comes from steam\n", "methane reforming (SMR) or electrolysis. It needs neither ecoinvent nor a\n", "commercial solver: the linear programs go to HiGHS and the cone to Clarabel,\n", "both installed with PULPO.\n", "\n", "As of now, PULPO considers uncertainty in the biosphere flows and the\n", "characterization factors only; technosphere exchanges are deterministic.\n", "\n", "1. the deterministic optimum, solved in reduced space (`formulation='reduced'`);\n", "2. the declared uncertainty and its closed-form moments;\n", "3. a chance-constrained front over the reliability level $\\lambda$, with an\n", " uncertain capacity of the electrolyzer;\n", "4. the exact variance decomposition at a solution;\n", "5. a screening of the parameters that declare no uncertainty;\n", "6. an out-of-sample validation of every point, and the check that the coverage\n", " is exact when the impact is exactly normal;\n", "7. the building blocks for formulations of your own.\n", "\n", "The method is described in the documentation (*Uncertainty*). The one plot\n", "needs `matplotlib`, which PULPO does not install; without it the plot is\n", "skipped and everything else runs." ] }, { "cell_type": "code", "execution_count": 1, "id": "7bc55a4b", "metadata": {}, "outputs": [], "source": [ "import contextlib, io\n", "import numpy as np\n", "import pandas as pd\n", "\n", "from pulpo import pulpo\n", "from pulpo.datasets.soc_demo_database import setup_soc_demo_db, METHOD_KEY\n", "from pulpo.utils import uncertainty as unc\n", "from pulpo.utils.utils import is_bw25\n", "\n", "pd.set_option('display.precision', 4)\n", "pd.set_option('display.width', 120)\n", "\n", "quiet = io.StringIO()\n", "with contextlib.redirect_stdout(quiet), contextlib.redirect_stderr(quiet): # it reports every write\n", " setup_soc_demo_db()" ] }, { "cell_type": "markdown", "id": "82ffc31c", "metadata": {}, "source": [ "## 1. The system and its deterministic optimum\n", "\n", "One unit of ammonia is demanded. Its hydrogen comes from SMR, which emits CO₂ and\n", "draws natural gas (with its methane leaks), or from electrolysis, which emits\n", "nothing directly but draws much electricity. The LCIA method is GWP100 over CO₂,\n", "CH₄ and N₂O." ] }, { "cell_type": "code", "execution_count": 2, "id": "0b42fb6c", "metadata": {}, "outputs": [], "source": [ "project = 'soc_demo_project_bw25' if is_bw25() else 'soc_demo_project'\n", "databases = ['soc_demo_background_db', 'soc_demo_foreground_db']\n", "method = str(METHOD_KEY)\n", "\n", "worker = pulpo.PulpoOptimizer(project, databases, {method: 1})\n", "with contextlib.redirect_stdout(io.StringIO()):\n", " worker.get_lci_data()\n", "ammonia = worker.retrieve_processes(processes=['ammonia synthesis'])[0]\n", "smr = worker.retrieve_processes(processes=['hydrogen SMR'])[0]\n", "electrolysis = worker.retrieve_processes(processes=['hydrogen electrolysis'])[0]\n", "\n", "choices = {'hydrogen': {smr: float('inf'), electrolysis: float('inf')}}\n", "with contextlib.redirect_stdout(io.StringIO()):\n", " worker.instantiate(choices=choices, demand={ammonia: 1})\n", "\n", "lci = worker.lci_data\n", "process_name = lambda j: lci['process_map_metadata'][j].split(' | ')[0]\n", "flow_name = lambda e: lci['intervention_map_metadata'][e].split(' | ')[0]\n", "SMR, ELY = (lci['process_map'][a.key] for a in (smr, electrolysis))" ] }, { "cell_type": "markdown", "id": "4ce50a22", "metadata": {}, "source": [ "With the technosphere matrix $A$ square and invertible, every scaling vector that\n", "meets the balances is $s = s_0 + S v$ with one free variable $v_k$ per alternative.\n", "`solve(formulation='reduced')` solves the LP over $v$ (here: two columns instead of six\n", "processes) and recovers $s$ with one forward solve. The full Pyomo model gives the\n", "same optimum." ] }, { "cell_type": "code", "execution_count": 3, "id": "85135307", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "impact, reduced: 1.6505600737 full: 1.6505600737\n", "SMR 0.0000, electrolysis 0.1780\n" ] } ], "source": [ "with contextlib.redirect_stdout(io.StringIO()):\n", " worker.solve(formulation='reduced')\n", "s_det = unc.current_scaling_vector(worker.instance)\n", "reduced_impact = worker.instance.impacts[method].value\n", "\n", "with contextlib.redirect_stdout(io.StringIO()):\n", " worker.solve() # the full LP, HiGHS\n", "print(f\"impact, reduced: {reduced_impact:.10f} full: {worker.instance.impacts[method].value:.10f}\")\n", "print(f\"SMR {s_det[SMR]:.4f}, electrolysis {s_det[ELY]:.4f}\")" ] }, { "cell_type": "markdown", "id": "d2ea0a65", "metadata": {}, "source": [ "## 2. Declared uncertainty and closed-form moments\n", "\n", "`import_declared` takes every nonzero entry of $B$ on a characterized flow and every\n", "CF of the method, with the distribution the database declares. Parameters without\n", "one are *undeclared*: they stay deterministic (mean = amount, variance 0) and are\n", "listed for the screening in section 5. CO₂'s CF is exact (a normal with zero\n", "width): CO₂ is the reference gas of GWP100." ] }, { "cell_type": "code", "execution_count": 4, "id": "85804ae9", "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", " \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", "
groupsubgroupstatusuncertainty_typen
0Ifsoc_demo_background_dbundeclared00
1Ifsoc_demo_background_dbdeclared23
2Ifsoc_demo_background_dbdeclared32
3Ifsoc_demo_background_dbdeclared51
4Ifsoc_demo_foreground_dbundeclared01
5Ifsoc_demo_foreground_dbdeclared41
6Cf('soc demo', 'climate change')undeclared01
7Cf('soc demo', 'climate change')declared21
8Cf('soc demo', 'climate change')declared31
\n", "
" ], "text/plain": [ " group subgroup status uncertainty_type n\n", "0 If soc_demo_background_db undeclared 0 0\n", "1 If soc_demo_background_db declared 2 3\n", "2 If soc_demo_background_db declared 3 2\n", "3 If soc_demo_background_db declared 5 1\n", "4 If soc_demo_foreground_db undeclared 0 1\n", "5 If soc_demo_foreground_db declared 4 1\n", "6 Cf ('soc demo', 'climate change') undeclared 0 1\n", "7 Cf ('soc demo', 'climate change') declared 2 1\n", "8 Cf ('soc demo', 'climate change') declared 3 1" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "data = unc.import_declared(worker)\n", "pd.DataFrame(unc.counts(data))" ] }, { "cell_type": "markdown", "id": "c959a104", "metadata": {}, "source": [ "With independent inputs, $X(s)=\\sum_e q_e \\sum_j b_{ej} s_j$ has, for every\n", "supported family (lognormal, normal, uniform, triangular),\n", "\n", "$$\\mathrm{E}[X] = \\mu^\\top s, \\qquad\n", "\\mathrm{Var}\\,X = \\sum_j d_j s_j^2 + \\sum_e w_e y_e^2,$$\n", "\n", "with $\\mu_j=\\sum_e \\mathrm{E}[q_e]\\mathrm{E}[b_{ej}]$,\n", "$d_j=\\sum_e(\\mathrm{E}[q_e]^2+\\mathrm{Var}\\,q_e)\\mathrm{Var}\\,b_{ej}$,\n", "$w_e=\\mathrm{Var}\\,q_e$ and $y_e=\\sum_j \\mathrm{E}[b_{ej}] s_j$. The second sum is the\n", "covariance that all processes emitting flow $e$ share through its one factor." ] }, { "cell_type": "code", "execution_count": 5, "id": "db484833", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "{'processes': 6, 'processes_with_variance': 4, 'uncertain_cfs': 1}\n", "at the deterministic optimum: E[X] = 1.8659, sigma = 0.8580, sigma without the shared-CF covariance = 0.8580\n" ] } ], "source": [ "mom = unc.compute_moments(data, worker)\n", "print(mom.summary())\n", "print(f\"at the deterministic optimum: E[X] = {mom.mean(s_det):.4f}, sigma = {mom.std(s_det):.4f}, \"\n", " f\"sigma without the shared-CF covariance = {mom.std_independent(s_det):.4f}\")" ] }, { "cell_type": "markdown", "id": "2118ddd1", "metadata": {}, "source": [ "## 3. Chance-constrained front\n", "\n", "The electrolyzer's capacity is uncertain: a triangular distribution between\n", "0.005 and 0.035 with its mode at 0.03. At a joint reliability $\\lambda$ the impact\n", "target and the capacity must hold together. Boole's inequality splits the failure\n", "budget $1-\\lambda$ across the $K=2$ events, which needs no assumption on how they\n", "depend on each other:\n", "\n", "$$\\min_s\\ \\mu^\\top s + \\kappa\\,\\sigma(s)\\quad\\text{s.t. PULPO's constraints},\\quad\n", "s_{\\text{ely}} \\le F_U^{-1}(\\varepsilon_1),\\qquad\n", "\\kappa = \\Phi^{-1}(\\lambda_z),\\ \\lambda_z = 1-\\varepsilon_0 .$$\n", "\n", "The capacity enters at its exact quantile. The impact row assumes $X$ normal\n", "(the only normality assumption; section 6 tests it). In reduced space the\n", "problem is a second-order cone with one column per alternative, solved by\n", "Clarabel." ] }, { "cell_type": "code", "execution_count": 6, "id": "0111d33b", "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", " \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", "
lambda_impactkappameansigmaadjustedsecondselectrolysiscapacity bound
lambda
0.500.7500.67452.15190.27402.33670.01190.01870.0187
0.700.8501.03642.15740.27262.43990.01550.01560.0156
0.900.9501.64492.16550.27202.61280.01570.01110.0111
0.950.9751.96002.16870.27222.70230.01480.00930.0093
0.990.9952.57582.17300.27302.87620.01500.00690.0069
\n", "
" ], "text/plain": [ " lambda_impact kappa mean sigma adjusted seconds electrolysis capacity bound\n", "lambda \n", "0.50 0.750 0.6745 2.1519 0.2740 2.3367 0.0119 0.0187 0.0187\n", "0.70 0.850 1.0364 2.1574 0.2726 2.4399 0.0155 0.0156 0.0156\n", "0.90 0.950 1.6449 2.1655 0.2720 2.6128 0.0157 0.0111 0.0111\n", "0.95 0.975 1.9600 2.1687 0.2722 2.7023 0.0148 0.0093 0.0093\n", "0.99 0.995 2.5758 2.1730 0.2730 2.8762 0.0150 0.0069 0.0069" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "capacity = {'uncertainty_type': 5, 'minimum': 0.005, 'loc': 0.03, 'maximum': 0.035}\n", "problem = unc.ChanceConstrained(worker, mom, upper_bounds={electrolysis: capacity})\n", "front = problem.solve([0.5, 0.7, 0.9, 0.95, 0.99])\n", "\n", "table = front.table()\n", "table['electrolysis'] = [p.s[ELY] for p in front.values()]\n", "table['capacity bound'] = [p.bounds[('upper', ELY)] for p in front.values()]\n", "table" ] }, { "cell_type": "markdown", "id": "88774a99", "metadata": {}, "source": [ "Diagnostics at every point: $\\sigma$ without the shared-CF covariance and the\n", "ratio $\\sigma/\\sigma_{\\text{indep}}$, the CFs' share of the variance." ] }, { "cell_type": "code", "execution_count": 7, "id": "82d748f8", "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", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
sigmasigma_indepsigma_over_indepcf_sharen_variablesn_rows
lambda
0.500.27400.27401.00020.001823
0.700.27260.27251.00020.001923
0.900.27200.27191.00020.002023
0.950.27220.27221.00020.002123
0.990.27300.27301.00020.002123
\n", "
" ], "text/plain": [ " sigma sigma_indep sigma_over_indep cf_share n_variables n_rows\n", "lambda \n", "0.50 0.2740 0.2740 1.0002 0.0018 2 3\n", "0.70 0.2726 0.2725 1.0002 0.0019 2 3\n", "0.90 0.2720 0.2719 1.0002 0.0020 2 3\n", "0.95 0.2722 0.2722 1.0002 0.0021 2 3\n", "0.99 0.2730 0.2730 1.0002 0.0021 2 3" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "unc.diagnostics(front, mom)[['sigma', 'sigma_indep', 'sigma_over_indep', 'cf_share', 'n_variables', 'n_rows']]" ] }, { "cell_type": "markdown", "id": "072003d3", "metadata": {}, "source": [ "## 4. Where the variance comes from\n", "\n", "$X$ is a sum of products of independent inputs, so its Sobol' decomposition stops\n", "at second order and is exact:\n", "\n", "$$S_1(b_{ej}) = \\frac{\\mathrm{E}[q_e]^2\\,\\mathrm{Var}(b_{ej})\\,s_j^2}{V},\\quad\n", "S_1(q_e) = \\frac{w_e y_e^2}{V},\\quad\n", "S_2(q_e,b_{ej}) = \\frac{w_e\\,\\mathrm{Var}(b_{ej})\\,s_j^2}{V},$$\n", "\n", "and $\\sum S_1+\\sum S_2=1$. No sampling, no confidence intervals." ] }, { "cell_type": "code", "execution_count": 8, "id": "989878be", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "sum S1 + sum S2 = 1.000000000000000\n" ] }, { "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", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
familyS1ST
parameter
Carbon dioxide, fossil @ hydrogen SMRsoc_demo_background_db0.79080.7908
Dinitrogen monoxide, fossil @ ammonia synthesissoc_demo_foreground_db0.12090.1209
Carbon dioxide, fossil @ electricity supplysoc_demo_background_db0.07160.0716
Methane, fossil @ natural gas extractionsoc_demo_background_db0.01440.0145
CF Methane, fossilCf0.00190.0020
Methane, fossil @ hydrogen SMRsoc_demo_background_db0.00010.0001
\n", "
" ], "text/plain": [ " family S1 ST\n", "parameter \n", "Carbon dioxide, fossil @ hydrogen SMR soc_demo_background_db 0.7908 0.7908\n", "Dinitrogen monoxide, fossil @ ammonia synthesis soc_demo_foreground_db 0.1209 0.1209\n", "Carbon dioxide, fossil @ electricity supply soc_demo_background_db 0.0716 0.0716\n", "Methane, fossil @ natural gas extraction soc_demo_background_db 0.0144 0.0145\n", "CF Methane, fossil Cf 0.0019 0.0020\n", "Methane, fossil @ hydrogen SMR soc_demo_background_db 0.0001 0.0001" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "point = front[0.9]\n", "dec = unc.decompose(point, mom, unc.families_by_database(worker))\n", "\n", "def named(frame):\n", " label = [f\"CF {flow_name(e)}\" if g == 'Cf' else f\"{flow_name(e)} @ {process_name(j)}\"\n", " for g, e, j in zip(frame['group'], frame['flow'], frame['process'])]\n", " return frame.assign(parameter=label).set_index('parameter').drop(columns=['group', 'flow', 'process'])\n", "\n", "print(f\"sum S1 + sum S2 = {dec.parameters['S1'].sum() + dec.interactions['S2'].sum():.15f}\")\n", "named(dec.top(6))" ] }, { "cell_type": "code", "execution_count": 9, "id": "72d35528", "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", "
S1ST
family
soc_demo_background_db0.87710.8772
soc_demo_foreground_db0.12090.1209
Cf0.00190.0020
\n", "
" ], "text/plain": [ " S1 ST\n", "family \n", "soc_demo_background_db 0.8771 0.8772\n", "soc_demo_foreground_db 0.1209 0.1209\n", "Cf 0.0019 0.0020" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "dec.families()" ] }, { "cell_type": "markdown", "id": "7d8437d2", "metadata": {}, "source": [ "## 5. Parameters that declare no uncertainty\n", "\n", "Undeclared parameters are deterministic in every result above. The screening\n", "asks what they would carry if they were uncertain: each gets the same coefficient\n", "of variation $r$, with its mean unchanged (a lognormal, mirrored for negative\n", "amounts), CO₂'s CF excepted. The ranking orders them by their deterministic\n", "contribution $\\mathrm{E}[q_e]\\,b_{ej}\\,s_j$, which fixes the order of their\n", "total-order indices up to the CF factor and does not depend on $r$." ] }, { "cell_type": "code", "execution_count": 10, "id": "387c7abe", "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", "
sigmaratioundeclared_share
r
0.10.27341.00520.0102
0.30.28441.04560.0844
\n", "
" ], "text/plain": [ " sigma ratio undeclared_share\n", "r \n", "0.1 0.2734 1.0052 0.0102\n", "0.3 0.2844 1.0456 0.0844" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "screening = unc.screen_undeclared(point, data, worker, exact_cfs=unc.co2_flows(worker), widths=(0.1, 0.3))\n", "screening.sigma" ] }, { "cell_type": "code", "execution_count": 11, "id": "25fec7bf", "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", "
amountcontributionST r=0.1ST r=0.3
CF Dinitrogen monoxide, fossil273.000.25941.0199e-020.0848
Carbon dioxide, fossil @ ammonia synthesis0.020.02005.3523e-050.0004
\n", "
" ], "text/plain": [ " amount contribution ST r=0.1 ST r=0.3\n", "CF Dinitrogen monoxide, fossil 273.00 0.2594 1.0199e-02 0.0848\n", "Carbon dioxide, fossil @ ammonia synthesis 0.02 0.0200 5.3523e-05 0.0004" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ranking = screening.ranking.copy()\n", "ranking.index = [f\"CF {flow_name(e)}\" if g == 'Cf' else f\"{flow_name(e)} @ {process_name(j)}\"\n", " for g, e, j in zip(ranking['group'], ranking['flow'], ranking['process'])]\n", "ranking[['amount', 'contribution', 'ST r=0.1', 'ST r=0.3']]" ] }, { "cell_type": "markdown", "id": "1b38e206", "metadata": {}, "source": [ "Width sensitivity: separate widths for each database and for the CFs. At the fixed\n", "decision the mean does not move, so the adjusted impact rises by\n", "$\\kappa\\,\\Delta\\sigma$; re-optimizing can only do better, which makes it an upper\n", "bound on the rise of the optimum. A re-solved point confirms it." ] }, { "cell_type": "code", "execution_count": 12, "id": "5342d053", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "re-solved at r = 0.3: the optimum rises by 0.02039 (bound 0.02039)\n" ] }, { "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", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
r[soc_demo_background_db]r[soc_demo_foreground_db]r[Cf]sigmaratiodelta_sigmabound
00.30.00.00.27201.00000.0000e+000.0000
10.00.30.00.27201.00026.6176e-050.0001
20.00.00.30.28431.04531.2331e-020.0203
30.30.30.30.28441.04561.2395e-020.0204
\n", "
" ], "text/plain": [ " r[soc_demo_background_db] r[soc_demo_foreground_db] r[Cf] sigma ratio delta_sigma bound\n", "0 0.3 0.0 0.0 0.2720 1.0000 0.0000e+00 0.0000\n", "1 0.0 0.3 0.0 0.2720 1.0002 6.6176e-05 0.0001\n", "2 0.0 0.0 0.3 0.2843 1.0453 1.2331e-02 0.0203\n", "3 0.3 0.3 0.3 0.2844 1.0456 1.2395e-02 0.0204" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "widths = [{'soc_demo_background_db': 0.3}, {'soc_demo_foreground_db': 0.3}, {'Cf': 0.3}, 0.3]\n", "sensitivity = unc.width_sensitivity(point, data, worker, widths, exact_cfs=unc.co2_flows(worker))\n", "\n", "wide = unc.compute_moments(unc.widen(data, 0.3, exact_cfs=unc.co2_flows(worker)), worker)\n", "resolved = unc.ChanceConstrained(worker, wide, upper_bounds={electrolysis: capacity}).solve_point(0.9)\n", "print(f\"re-solved at r = 0.3: the optimum rises by {resolved.adjusted - point.adjusted:.5f} \"\n", " f\"(bound {sensitivity['bound'].iloc[-1]:.5f})\")\n", "sensitivity" ] }, { "cell_type": "markdown", "id": "43d475a4", "metadata": {}, "source": [ "## 6. Out-of-sample validation\n", "\n", "Each solution is held fixed and every declared input is drawn from its own family\n", "(not a fitted normal), together with the capacity: $N = 200{,}000$ draws, each\n", "parameter from its own stream keyed by one recorded seed, the same draws for\n", "every decision. Per point: the coverage of the\n", "impact row against $\\lambda_z$, of the capacity against $1-\\varepsilon_1$, and of both\n", "together against $\\lambda$, with Wilson intervals; the normal-approximation error\n", "of the $\\lambda_z$ quantile; skewness; the empirical and the Gaussian CVaR; and the\n", "regret against the lowest level and against the best point of the front in each\n", "draw. The deterministic optimum is priced on the same draws." ] }, { "cell_type": "code", "execution_count": 13, "id": "e389f68b", "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", " \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", " \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", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
lambda_impactcoveragecoverage_lowcoverage_highnominal upper:4coverage upper:4joint_coveragequantile_errorskew
designlambda
front0.500.7500.76180.75990.76360.7500.75070.5705-0.01110.3829
0.700.8500.85130.84980.85290.8500.84950.7229-0.00180.3661
0.900.9500.94160.94050.94260.9500.94990.89440.02550.3542
0.950.9750.96600.96520.96680.9750.97530.94220.04490.3529
0.990.9950.98970.98930.99020.9950.99480.98460.09460.3535
deterministic0.500.7500.77550.77360.77730.7500.00000.0000-0.08421.7324
0.700.8500.80360.80180.80530.8500.00000.00000.21001.7324
0.900.9500.84250.84090.84410.9500.00000.00000.88491.7324
0.950.9750.86020.85870.86180.9750.00000.00001.35201.7324
0.990.9950.88870.88730.89010.9950.00000.00002.53191.7324
\n", "
" ], "text/plain": [ " lambda_impact coverage coverage_low coverage_high nominal upper:4 coverage upper:4 \\\n", "design lambda \n", "front 0.50 0.750 0.7618 0.7599 0.7636 0.750 0.7507 \n", " 0.70 0.850 0.8513 0.8498 0.8529 0.850 0.8495 \n", " 0.90 0.950 0.9416 0.9405 0.9426 0.950 0.9499 \n", " 0.95 0.975 0.9660 0.9652 0.9668 0.975 0.9753 \n", " 0.99 0.995 0.9897 0.9893 0.9902 0.995 0.9948 \n", "deterministic 0.50 0.750 0.7755 0.7736 0.7773 0.750 0.0000 \n", " 0.70 0.850 0.8036 0.8018 0.8053 0.850 0.0000 \n", " 0.90 0.950 0.8425 0.8409 0.8441 0.950 0.0000 \n", " 0.95 0.975 0.8602 0.8587 0.8618 0.975 0.0000 \n", " 0.99 0.995 0.8887 0.8873 0.8901 0.995 0.0000 \n", "\n", " joint_coverage quantile_error skew \n", "design lambda \n", "front 0.50 0.5705 -0.0111 0.3829 \n", " 0.70 0.7229 -0.0018 0.3661 \n", " 0.90 0.8944 0.0255 0.3542 \n", " 0.95 0.9422 0.0449 0.3529 \n", " 0.99 0.9846 0.0946 0.3535 \n", "deterministic 0.50 0.0000 -0.0842 1.7324 \n", " 0.70 0.0000 0.2100 1.7324 \n", " 0.90 0.0000 0.8849 1.7324 \n", " 0.95 0.0000 1.3520 1.7324 \n", " 0.99 0.0000 2.5319 1.7324 " ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "validation = unc.validate(front, problem, data, n=200_000, seed=2026, designs={'deterministic': s_det})\n", "columns = ['lambda_impact', 'coverage', 'coverage_low', 'coverage_high', f'nominal upper:{ELY}',\n", " f'coverage upper:{ELY}', 'joint_coverage', 'quantile_error', 'skew']\n", "validation.table[columns]" ] }, { "cell_type": "code", "execution_count": 14, "id": "fb33a13e", "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", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
thresholdcvar_gaussiancvar_empiricalregret_vs_lowest_lambda_meanregret_vs_front_best_mean
lambda
0.502.33672.50022.51330.00000.0143
0.702.43992.58112.60650.00550.0198
0.902.61282.72652.78170.01350.0278
0.952.70232.80512.88060.01670.0310
0.992.87622.96263.09040.02100.0352
\n", "
" ], "text/plain": [ " threshold cvar_gaussian cvar_empirical regret_vs_lowest_lambda_mean regret_vs_front_best_mean\n", "lambda \n", "0.50 2.3367 2.5002 2.5133 0.0000 0.0143\n", "0.70 2.4399 2.5811 2.6065 0.0055 0.0198\n", "0.90 2.6128 2.7265 2.7817 0.0135 0.0278\n", "0.95 2.7023 2.8051 2.8806 0.0167 0.0310\n", "0.99 2.8762 2.9626 3.0904 0.0210 0.0352" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "validation.table[['threshold', 'cvar_gaussian', 'cvar_empirical', 'regret_vs_lowest_lambda_mean',\n", " 'regret_vs_front_best_mean']].loc['front']" ] }, { "cell_type": "markdown", "id": "f5500b1b", "metadata": {}, "source": [ "The impact is right-skewed here (the lognormal CO₂ emission of SMR dominates its\n", "variance), so the Gaussian quantile sits below the empirical one at high\n", "reliability: the impact row holds slightly less often than it claims, the joint\n", "event falls short of $\\lambda$ by a similar margin, and the empirical CVaR exceeds\n", "the Gaussian one. The capacity holds at its nominal level, since it enters at its\n", "exact quantile. The deterministic optimum was solved without the capacity and\n", "runs the electrolyzer at 0.178, above the capacity's maximum, so the bound never\n", "holds." ] }, { "cell_type": "code", "execution_count": 15, "id": "f32e836f", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAk0AAAFUCAYAAAA5ywiWAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjcsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvTLEjVAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAYWtJREFUeJzt3Qd4VMXXBvA3vVeSEEgCoRfpvQooRUGwoYgFBMSChfK3oQgqKiqIICJ8gCgWBFGaiBRBepAmgvTeE9J7z37PmWWXJCRh07a+v+e55G6fvTdkz545M2On0Wg0ICIiIqIS2Zd8MxERERExaCIiIiIyEDNNRERERAZg0ERERERkAAZNRERERAZg0ERERERkAAZNRERERAZg0ERERERkAAZNRERERAZg0ERERERUGUHTtm3b0L9/f1SvXh12dnZYuXLlbR+TmZmJt99+GzVr1oSLiwvCw8OxcOHC0r40ERERkck4lvYBqampaN68OYYPH46HHnrIoMc8+uijiIqKwtdff426devi2rVryMvLK0t7iYiIiCwjaLr33nvVZqh169Zh69atOHv2LPz9/dV1kmkiIiIisuqgqbRWr16NNm3a4NNPP8X3338PDw8PDBgwAJMnT4abm1ux3Xmy6UhWKi4uDlWqVFFdgkREREQVRaPRIDk5WZUe2dvbmy5okgzTjh074OrqihUrViAmJgajRo1CbGwsvvnmmyIfM2XKFLz33nuV3TQiIiIivUuXLiE0NBTFsdNIeFVGkvWRQOiBBx4o9j69e/fG9u3bERkZCR8fH3Xd8uXLMXDgQFUfVVS2qXCmKTExETVq1FBvxtvbu6zNJSIiI5O/8/LtXVy9elX1NhCZm6SkJISFhSEhIUEfq5gk01StWjWEhIQUaESjRo1UKuzy5cuoV6/eLY+REXayFSYBE4MmIiLL4eDgoN+Xv98Mmsic3a4EqNLnaercubP6dpGSkqK/7uTJk6rPsKQUGBERWT4JkuRLsmwMmMjSlTpokuDn4MGDahPnzp1T+xcvXlSXx48fjyFDhujv//jjj6sC7mHDhuHo0aNqnqfXXntNTVlQXCE4ERERkcUHTfv27UPLli3VJsaNG6f2J06cqC7LHEy6AEp4enpi48aNqp9QRtE98cQTanLML774oiLfBxEREVGlKlchuDELtKQmSgrCWdNERGWVm5uL7OxsHkAjkkE9r7/+utqXqWeKqlclqmxOTk4F6uvKGmdUeiE4EZGpyXdDGcErGW8yLpln78EHH1T7MvinpDlwiCqTr68vgoODyzXfI4MmIrJ6uoApKCgI7u7unCTXyNm99PR0/WoQJX3bJ6qsL01paWm4fv26flR/WTFoIiKr/9DWBUwyKIWMf/x1ZJJjBk1kCrqBZxI4yd+Csv4eMk9KRFZNV8MkGSYisl3uN/4GlKeukUETEdkErltJZNvsKmDtWgZNRERERAZg0EREJjV++WG1kWG6d++OMWPG8HARmQALwYnIJhk7UJvyUFPYIpliQDcZMqcbIEvHTBMRkZXIysqCOdaRyEgl2VhXRpaOQRMRVWq3G7veyi41NVWt5SnLUcncMp999lmB22Xeo8mTJ6v7yCzGzz77rLr+jTfeQP369dVoodq1a+Odd97RjxiSGY8lgJElsXSTT/r7+6NDhw765/3hhx8QFhamD8Reeukl9foyZUDNmjUxZcqUcrwrIsvFoImIKh2Dp7KRxc23bt2KVatWYcOGDdiyZQsOHDhQ4D7Tpk1D8+bN8c8//6jgSHh5eeHbb79Vi6TPnDkT8+fPx+eff65uk6UiWrRooZ5LHD58WGWA5PGyILuQ1+zWrZval3VCV69ejZ9//hknTpzAjz/+qII1Q0lQJgu7yyb7RJaMNU1EZDTMOhlOApivv/5aZX3uvvtudd2iRYsQGhpa4H533XUX/ve//xW4bsKECfp9CXBeffVVLFmyRL8GnBSTS9Ak18vPXr164fjx49ixYwfuuecedZ3uvrIAe7169dClSxcVXEmmqbSzMcfGxqr9GjVqlOqxROaGQRMRmV1AZatF0/mdOXNGdY21b99ef510ozVo0KDA/dq0aXPLY5cuXaoyRPIcEnzl5OQUWIRUskgSkMls3ZJV6t27t1qTS4KlZs2a4fTp0yqwEk8//bQKquR1JaC677771P2JbBG754jIbLEu6vY8PDwKXI6IiMATTzyBvn37Ys2aNarb7e233y5QJH7nnXciOTlZdfVt27ZNBUi67JMEUdWrV1fZJdGqVSvVtSa1U7KG3KOPPoqBAwdW+LkmsgTMNBGR5dBogNQYIP48kJUsY7MAB2fALxzwri5DtWAt6tSpAycnJ/z999/6bq34+HicPHlSX29UlF27dqkuNAmUdC5cuHDLau+SUfryyy/VazRs2FCtxzVo0CAVaBV+fslSyW2yScAkGae4uDiV+SKyJQyaiMis+WRFolHSTmDxh8CFnUBmUtF3dPIAQloB9XoB9e8FAuvDksmIuREjRqhicFloWIIaCYRuN9eRZIikDklqmNq2bYvff/8dK1asuOV+klmaNWuWPmskAVCjRo1U197s2bP195s+fboaOSdzLclrL1u2THXlSeBFZGsYNBGReWaUzm7Bk+c/UwGTPfKPurIDvEMAN1/t/bJTgYRL2p/nt2u3jROBGp2ADs8DNbVF1JZo6tSpqiapf//+akScFHzLlAElGTBgAMaOHaumCcjMzES/fv3UqLp33323wP0kmzRjxgx97ZKQ/X///bfAdfK6n376KU6dOqWmKpBAbO3atZyokmySnUaGNpi5pKQkNUxW/ljkL2YkIusbKReSdhz9rs1CrdR/9ded82iOE14dcMqzHV4edB/g5FrgNabc3xCIPQOc2wacWq8CLuTlqNsyQu/EuXbvo1b9RmqeITIuKTaXuioh2SoJvIhMISMjQ9Xn1apV65a/BYbGGcw0EZFZBEsuuanod3UW2savUZez7Fywz/8+7K7yIKJd880LdCNgKsDBCQhqqN3aPwskXQX2fg3s+xpIOAekXgfiPXDNJRi5do4I9XMv8/uj0pEuPZlHSrdPZMkYNBGRyYWn/otHLn4A/+xr6vIB3z7YEPwsEp2rlu0JpSj87neATi8DO+dou/SyUxCkOYcE52AADJqMReZ2kmJzImvAoImIykzfPVbWeZU0GnSN+Qn3XJsDe2gQ51QNy2pMwHkPbWaiIto2qe9Y4PRJwNERDshElayrQEIW4BMC2DHzQUSGY9BERCbpknPIy8IDV6ahTfxadfmA3z1YXX0sMh0KzjtUIaT7zq8mkpJj4JUTC7u0GCAnE/APB+z5Z7AyydIply5dUvuynh276MiSlfprlkyEJiM5ZPIzSbuuXLnS4Mfu3LkTjo6Oat0jIrJdzrlpGHbufypgyoM9Vlcfg2WhbxsUMJV2wsvpG04iPi0LkUmZSHIKQKzzjQyTzPMUcxLIuTnpI1U8GWsUHR2tNgsYd0RUsUGTrLotRX355/EwREJCglqJW7eGEhHZJtfcZAw/Nw51Uv9Bhr07vq01DREBA402MWWGgydQpR5g76TNNsWe0v4kIrqNUuel7733XrWV1vPPP4/HH39cDTctTXaKiKxn0V3XnCSMODcWoeknkObghW9qTcdl90bGb5uzOxBQD4g9DeRmaX9WqQs4ulRYW4jI+hilCvKbb77B2bNnMWnSJGO8HBGZIae8DAw9/4YKmFIcfDG/9hcVGjCVmgRIEjg5uNwMnHKzTdceIjJ7lV4BKbPIvvnmm9i+fbuqZzKEzGIrW/5Jp4jIctlrcjD4wkSEpx1Gur0nvq49A5FudU3dLO26dRI4SW2TCpzOAAF1WRxORMbPNMlMsNIl995776F+fcPXgZoyZYqamVO3yYgLIrJQGg0euDwNjZJ3IdvOGYtqfVLhAVNpCsOLHFknXXMyii4nHYg7B2jyL9tChgoPD1dLs1SUp59+Gg888IBJT4AhbUhLS8PDDz+sZpKWAVJSw0vWqVKDpuTkZOzbt0+tgSRZJtnef/99tbaR7G/evLnIx40fP15NZa7bdMNVicjydI75Wc3yLaPkFtecjAsVMAdTpXTV+de5MaouBUi8rF3XjsgAixYtUr0pu3btwrVr19SXfapYW7ZsMYuAtFK75yTqPny44Le/r776SgVLv/zyi1r/pSguLi5qIyLTyp+9KcsElvWS/0bfa9qRtmurvYTj3p1htqQ43C8ciDsLpMUCTm6ARyCsSVZWFpydnY36mjIvU9Om2t8da52j6cyZM2jUqBGaNGlSrp4ZCQqs9RhZi1KfHVlx++DBg2oTsvid7F+8eFGfJZKpBdST29urX6L8W1BQkFooT/Y9PCphEjsiMgv+mVcw+MIk2CMPe/36YWfAIzB7rj6AV3XtvmSbMpNN1pTu3bvjlVdeweuvvw5/f38EBwfj3XffLXAf+bt7//33w9PTU31JffTRRxEVFaW/Xe4v8+ItWLCgwCKl8uH8f//3f7jvvvvg7u6uPvAjIiJw+vRp9bryt7lTp04qGNCRfXmtqlWrqtdr27Yt/vzzz9u+D3kt3Rdh2S/syJEjqh3Sfi8vL3Tt2rXA6+Ynta5yTHSfI126dMHevXv1t8fHx+OJJ55AYGAg3NzcUK9ePTUQSUd6LeQY+fr6qmMq7+f8+fMFApdx48ap26tUqaKO/e3mlpLj9dlnn6k5DOX9yWVdW+Sz0M/PTx1jGXUuNb463377rXqd1atXo3Hjxur4yPmULs4PPvhAPVaOc82aNdV9ZJ4r3blu1qyZ6sW53XEv7TkWc+bMQZ06dVRw3aBBA3z//fcV8ryrVq1Cq1at1HmrXbu2KtvJyckp8Lzye/rggw+q55VzJ+9byDnq0aOH2pfjKfeVblMhCRgJyuV8yznr2bOnmhrJbIImOVGyUrVsQn7BZH/ixInqsqQmdQEUEdkmx7xMPH7hHbjlpeCCexOsCvmf0eZhui35EMxKLX5z9tB212WnA9ePA+kJJd+/NFspu/yk20c+hP7++298+umnqrxh48aN+pm25UM0Li4OW7duVdfLKOVBgwYVeA75MPv111+xfPly/ZddMXnyZPXBLNc1bNhQ1Z8+99xz6ouv/J2XYEFKK/J/Ye7bty82bdqEf/75B/fcc4+a6Lg8f++vXLmCO++8UwUM0gOxf/9+DB8+vMCHaX4SxMh7keNy4MAB1K1bF3369FHHQLzzzjs4evQo/vjjDxw7dkwFAAEBAeq27OxsdV8JzKQrTSZblgBE3odk4IQEPxLMLFy4EDt27FDPu2LFihLfgxzXkSNHomPHjurzTy4L+VCX4ygf/BJUyPGU4yftyF8L9cknn6hgQYJHCQbF559/js6dO6vj3K9fPzz11FPqXD355JPqfUtQI5dvF9CV9hzLex09ejT+97//4b///lP3HTZsGP76669yPe/27dvV/eW55fxI0CXH+cMPPyzwvBJISVB76NAhdawkAJZzIHXNct7FiRMn1HGeOXOm+jl48GD1OyPnW7rwHnrooUqdRNVOYwFTtMroOekjlvom+TZCRKbrnjOk4Pr+K9PQIXalmlpgVr2FSHLWfhiYgpejBt1DNKgeWgOOzi6wy05DyFd1TNOYt65qgzIDyLd2yXzIB45Ou3btcNddd+Hjjz9WQZJkLyTbrxssIx9Id9xxB/bs2aMyQZJp+uijj1RwItkXHfmmPmHCBPXhJ3bv3q0+9L/++mv1ASSWLFmiPjDT09OLbaP0GMgcfLoPSMmSjBkzRm06EtzJ64uQkJAC3U9vvfWWeh35ICxqUV8JPKSGReb2k+yBZBnkw1Y+pIUEILrXfO211zBgwAAVJEnQU9gPP/ygMjjy4arLeEmwJNkeef7evXurlS7Gjh2rnktI8CYZutatW5c4v6C8vgQQ8qEtJKMkg58kMJOsi4iNjVXnSQK+Rx55RL0POb7yOJkwWkfej2TbdBmeyMhIVKtWTQWEEjTnP18SNEgGsihlOccSqMnvz7x58/TPI0GMHPvff/+9zM8r2R+Z2FqCqvznQ4Lgq1evFvm88poS1EoALIGtHFvJNkkGT86ZkABSzo1koiQjdzsZGRnq/0v+rGtp4wx2nhJRhWqasEkFTEIW3zVlwGTppBsmP/nwvH79utqXD3/5EM4/uli6eeQDRW7TkQ+T/AFTUc8tXW5CV3uku04+ZHRTvkim6dVXX1XdMfIa8oEmr3O7TJN8L5cuQ9kKf0eXgEEChKICpsKku0eCJPlg15HHSSCpe78vvPCC+sCWLkn5QJbCbB0ZgCRZN8k0Sdtlky46eY/y3PJhKUFI+/bt9Y+RAUtt2rTRX/7xxx/1j5Utf0Cbn7RHHpv/uaTrSLq78p8b6QIrfI4NPTdC97tQnNKeY2lb/uMr5HL+NpfleeXYS8CX/9hJdk6Ot2TbinpeybBK8FLSe5RgU4IxeW0JROfPn6+CqsrElSqJqMJ4Z13Hg1emqv2/Ap/CSa8OZnd0NY5u2oyPIdLigMQbo3dl6RUpFi8Pp9I9vnAwId/GJXNTGsXVjuZ/bl3mpajrdK8nAZNkt6ZNm6a6xaSGZODAgfqurbKQ56hIknm7cOEC1q5dq9oqH6gvvviiarMEfZKVkMCnsKKCyqJIJit/ICSZs/KQ919UnVdpz01xKup5yvu8KSkpqutNus4Ky5/xKe3vu6wwIudZguMNGzZg1qxZePvtt1V3dnEDzcqLmSYiqhB2mjwMvDwFbrkpuOzWEH8GjzDPIyt/0KWLrNB2OdVOv+mv9w3TFobLSLrUaMDRtcjHGrxVYF2XZHyksDn/lCzSPSfdWZJxqmjS1STdZVKoK9/spVsofxF1WUhmQbI1+et8iqMrTpZ26MjjpBA8//uVAGjo0KGq+0fmjNJ1NUkRsnSbSd2QBH35N92cgJLJkw9cHemekzorHclS5X9ccUGfnBt5bP7nku456YasjHNTUaTd+Y+vkMvlbXOrVq3Uey983GUzdLSgbtSndFkXDqwkGyZBmdSAyf1uV4dWHsw0EZFBblfL1CF2Oeql7EWWnQuWhk1Enp35/nm5HH+zSyDU7zbZH59Q7dxNuZlA8lXAxzwm25U6EQlepFhWggP5kB41ahS6detWoEuposhoJilyluJv+aCSGpvSZigKk1ooyQ489thjqt5FAhepkZEuN+nKKpwxk+43qTeSbrUaNWqo4njp3hkxQhugy4AkySZJXY6MtFuzZo0KBIQcp6lTp6rieekqCg0NVVkpeU/SlSeXpVBZ6sXkvUqB8/Tp08s0L5A8Xl5HuqCk6FmCLVkZQzJTcr25kmMrNUwyuEt+v3777Td1fAwZJVkSOS8y2k7OmWQnJVCSLjspNpc6M0NIN7P83sk5lSJxCVileF4GJkg9mgTDEqTKKEPdOa8MzDQRUblVybyMe699pfbXVRuFGNcaFnNUJYDKH0Tdcp2DI+B74/2kxph0GoL85ANEhnFLcbSMQJMPORnKvXTp0kp5PQkg5LWksFkCJxmJJhmEkkhmSgrXiyN1PjJqTrpvJNiTgEfqUoqrcZKARmbeltFk8tpSo7R+/XrVLiFZBgm+JIMlx0S6b6TGScgwdpkWQD64pZtIPlgl2JLaG13hr4wak+eWTJUUN0uwI5m1spCpDuT9SLAgzyX1XNJtaEj9lqnIzOcyKk26MyXwlIBP3oduGoWykt8VCXakC00GKHTo0EGNEDSkeFtHAk7JJknwKTVTEnDLeZNzKkGUFN5LIbmMgJRu2srC0XNEdItSLUmi0WDEuTGom7Ifpz1bY2Gtz6GRmbXNROHRc6VRIAuVcAlIi9GuVxfYELB3qPjGWhkJhGSTWiAhGQwJZIhMgaPniMjkWsX/oQImWVduZchrZhUwVSjv6oC9k3Zh35SbE0hS0WQ0moxKk+wNkbUw36IDIjJ7Hjnx6HftS7X/Z9XhiHUJhdWSzJLUM8Wf1QZNbn7aAnEqktQnXb58WXVLSVePOoRcIoQsnJV+JSQiY7jv6ky45ybhqmtd7Ah8zOoO+i31Tm4+2qVWRMJFLuprYO2VFO0WN7yeyJIwaCKiMqmfFIEWCX8iD/ZYHvqmWY+Wq1Ayms7OAchO0xaGE5HNYNBERGVaW27A1c/V/q6AR3DFvaHtHEUpBPeupt2XKQikxomKJdMSyFIZspV3igIiU2PQRESl1jnmZ1TJuopExwBsrGqmk1hWJvcA7ezemjwg6ZqpW2PWpKZJFzRZwFKnRCVi0EREpeKdHY0e179T++uqvYAsh3IuLWKJpDZHuulEehyQlWrqFhGRETBoIqJSuefaHLjkpeOCexMc9O1tu0dPlkVx89fuJ15mUTiRDWDQREQGq5F6GC0TNiAPdvit+pgKXUvNYuduknmppChcMk5EZNUYNBGRwQvy9r86Q+3v9+trW8XfxXFwAjyDtftJV4G8gouJEpF1YdBERAZpGb8OoeknkGHvgQ3Bz9nUUdPN11R4jTrFMxBwcAHycoCUSKO2S1dY/e677xa4TESVg0ETERk0xUCvqAVqf3PQUKQ43ajlIW33nE+I9kikRAM5xpuCYM6cOZg3bx5SU1PVQqayeCkRVR4bmY2OiMqzOG+nmGXwzb6OeKeqiAh4mAezMBdvwNkTyEoBkq8Bfoav3l4eo0aNwtSpU/HFF19g8+bN6Ny5s1HPjQRqslr9ww8/jMWLFxd5H1k6pVGjRvp9IkvG32AiKpF7TiK6X/9B7W8MfgY59i48YoVJQbwUhQspCM9ON8oxmjt3rlrj7ZVXXsFvv/2G7du3l/s5u3fvjjFjxhh03/Hjx+Ozzz7DTz/9hNOnTxd5H1k6xcPDQ21lXUZl9uzZCA8Ph6urK9q3b489e/aUeP/k5GT1HmrWrKmWb+nUqRP27t1bptcmyo9BE5GNZZV0m6G6Xf8ebnkpuOZax7anGDBkCgLdunRSFG4Ezz33HJ599lkVkHz88cfo0qULjEkCthEjRqgM0uHDhv9OlcbSpUsxbtw4TJo0CQcOHEDz5s3Rp08fXL9+vdjHPPPMM9i4cSO+//571a7evXujZ8+euHLlSqW0kWwHgyYiG2VI8OSbFYlOsb+q/XXBz0Mja65R8bxuZJsyk4DMlHIfqY8++khlZwpvM2ZoRzHqMje6QnBDMjnr1q1TwZWvry+qVKmC++67D2fOnFG3Pf3009i6dStmzpypf63z58+X+Hw5OTlwd3fHf//9V+TtsnRKZGSk2sqyjMr06dMxcuRIDBs2DI0bN1bZNXm9hQsXFnn/9PR0/Prrr/j0009x5513om7duur4yE+pASMqDwZNRFSsXpHz4ajJxhmPVjjp1YFH6nacXLVLrIikK+We8PLll1/GtWvX9JsED9LlNHDgwDI/pxSNS+Zm37592LRpk8oSPfjggyqgkWCpY8eO6nV0rxkWFlbi802YMAEpKSnFBk0yom/y5MmoXbs2vL294enpWeR28eLFWx6blZWF/fv3qyyRjrRXLkdERBQbxOXm5qquvPykm27Hjh0GHiWiCioEl9EZUngov8jyH2rFihV44IEHir3/8uXLVXR/8OBBZGZm4o477lBRv6RXich8Vc04ixYJG/TLpVjbRJZpqcUvfWLv4FDgQzf/fVOdCwZC8iEuH8j62+29gPQrQFo6YHcVcPNV10sXWml5eXmpTbzzzjvYsGEDtmzZgtDQG0u4lIEUbecnGZvAwEAcPXoUTZo0gbOzs8rkBAffmH+qBPI5IJmffv36FRs0iYceekgFOpIpcnAoOltZvfqNLF0+MTExKgCqWrVqgevl8vHjx4t8HjleEvhJoCYF6HJfqbmSIEuyTURGDZrkW4r0KQ8fPlz9RzAkyOrVq5dKM0s6+JtvvkH//v3x999/o2XLlmVtNxFVsrujFsIeGhz26Y7L7trRT9akQVhQsbfd1asPFi1drr/cokE40tOKmKMJQIfOXRGx4+ZQ//A69dSHfWHlmUNp4sSJqj5HAiYpiC6PU6dOqeeTv8HSTl2XmWR6JGgylDxOaqpeeuklVZz95JNPIjs7G05OTkXWPskmQUtxQVNFkmMln1EhISHq9Vq1aoXBgwerII/IqEHTvffeqzZD6fredSR4WrVqlRrpwaCJyDxVSz+Fpolb1HIpf1Ydburm2DQpgP7uu+8qJGAS8qVVuvjmz5+vsjsS/EiwJF1hpTFr1iwVdL3//vsq4JKASbI/TZs2veW+8mVZtpKmHJBMV40aNQpcFxAQoIKeqKioAtfL5ZIyYXXq1FG1WfIlPykpCdWqVcOgQYNUFyGRRc3TJP9BZTiov3/xk+NJN55sOvJLT0TG0zNSO5HlId+7cd3VOj9oTly6XmL3XH4HTxRfDG1XKBDQF04nRwIpUdrZwgMblDlgWrRoUYUFTLGxsThx4oQKmLp27aquK1znI91z0iVWEhmFJt2F0u0l3Y716tWDi4uL6qIrKmgqa/ectKV169aq9kpXBiKfIXJZMly3o5vqID4+HuvXr1fF4UQWFTRNmzZNFQ0++uijxd5nypQpeO+994zaLiLSCk07hsbJO5EHe2wKst4sk3spaoxKc1997ZJbTUCTCmhyAfvSzxL+wQcfqHrQ1atXq/oqGX0m/Pz8VIBSFvJYGTEns4hL9kUyRDJBZX4SnEnXnQR/UqAtX3ALZ4hkXijpcZBaJuHo6Kjqh4qraypP95wUrQ8dOhRt2rRBu3btVO+FZJBkNJ348ssvVW2tBFI6EiBJd2iDBg3U/FGvvfYaGjZsqH8MkUWMnpMZYyUY+vnnnxEUFFTihGmJiYn67dKlS8ZsJpFN63ljuZSDfr0R41qwu4RKwd4R8Ay6mXUqRU2TfODLgJvo6GhV1CwBjm47dOhQmU+DBD9LlixRtT3SJTd27Fj1Ovm9+uqrKrCRrJAUiBce1bZmzRo1+7iMtMtPMkwlFYOXlXSryZdtqcNq0aKFGlQk0yboisOli1A3ZYKOfG68+OKLKlAaMmSImmJBAqmi6q2ISsNOU47qRJnD43aj53TkP6oU5i1btkz/7cRQ0j0n31LkP4IMWSWisrndvEw1Ug/jhTMvIBcOmN5gMeJcbqypZsG8HDXoHqJB9dAacHSunNnMQ/3ci74hLxeIOqLNNvnWANyrwNbIR4z0LgjJXJV1VnCi8srIyMC5c+dQq1atW6akMDTOMEqmSfq9JS0qP0sbMBGR8fSM+lr9POB3r1UETCZn7wB4Vc2XbSr95I6WToIk3dQJDJjI5mqa5BtD/jWGJGqTdKn0e8vIB+lakyJBGe2h65KT/mhJ5cqwVF2/vMxrIlEdEZmHsLQjqJeyT2WZNlcdaurmWA+Z7DLlOpCbBaTH22S2ichalDrTJLPIylQBuukCpEhP9qW/WciEl/n7wKXgUGZolf7l/P3yo0ePrsj3QUTl1D1K+0XnH78+SHCuxuNZkdkmD11tU1S5Zwm3NDLaTdaJk60sy6gQWXSmSVbALqkM6ttvvy1wWYbKEpH5z8ukGzG3NehJUzfH+nhItikKyM0EMhIANz/YCvm80H2RlpF7RJaMa88REbpf12aZDvv0QIwLR8xVTrYp0GazTUTWwujzNBGReY2aC8y4gCaJ2ozwX0FDjNgqGyNBU+p1ICcdyEwCXFnTSWRpmGkissIA6XZTC+TXPfp7tcbcEe+uiHKrA2ujy+lo9HsV73J8mtpK5OCo7aYrw7xNRFR+5Vn/UYdBE5EN88u8gubxG9X+FivNMmXkArl5GuRkZlT6a+mCp2IDKFUQbgdkpwFZ2rmLiMg40m4sul2eSU7ZPUdkpQzJNnWLXgwH5OKkZztcdm8Ea5SjscPZRMDZIQay4qWjiyvsJHCpZBkZxXwndfLRTj0QdwXwK/96cuYu/zp2MrlgaZdRIaqIDJMETDKC09fXt1y/gwyaiGyUd3Y0WsevVft/Wfm8TMeTJIDJQ+3c63CwN86M1OlxzkXfkJcDJMVoOw5jMiSKgzWTaQZkqRMh69kVXseOyFgkYAoODi7XczBoIrJRXaN/gqMmG+c8muO8R3NYNzscT3LA6WQNXB1UB1mlG9e7VvE3bloCHFsF1OwM9C+4hpu1kXn6dOuH1qxZUy3uS2Rs0iVXEVlO/vYS2SC3nCS0jfvN5kbMSVddSo5xXqvw2lYFdBgO7PsKOLIE6PoyENwE1uyee+4xdROIKgTzpEQ2qH3sCrjkpeOqa12c8mxn6ubYnip1gMb3a/d3fWHq1hCRgRg0EdkYx7xMdIr9Re1vD3xcVlQ1dZNsU6dXtD//+xVIvAxrlZ2drVaKkE32iSwZgyYiG9Myfh28cuKR4BSEQ753mbo5tiukFRDeVVsYvnsOrFVWVhaGDRumNtknsmQMmohsiJ0mD12jl6j9HQGDkGfHskaT6nxj4fL93wLpCaZtCxHdFoMmIhvSKGkHArMuId3BE3v9+5u6OVS3JxDUWDvR5b6FPB5EZo5BE5ENuTN6sfq52/9BZDm4m7o5JPVknV7WHoe/5wI5mTwmRGaMuXkiG1Ez9RBqpv2HHDsn7AoYaOrm2LT8s7VPGTAQ2DQZSL4KHFoKtLKdKSCILA0zTUQ2lmU64HcPUpyqmLo5NreAcrHL2jg6Ax1e0O7vmiVTaBu1fURkOAZNRDYgMOMCGiftQB7ssD1gsKmbQ4W1fhpw8QZiTgKn1vP4EJkpds8R2YAuMdoRc8e8uyDGtYapm2MTisosFZttcvUG2gwDds7Ubg3uhbVwcXHBzz//rN8nsmQMmoisnFd2DFrFr1P722QySzJP7V8AIr4CLkYAl/YCYW1hDWStuUceecTUzSCqEAyaiKxAsRkMAB1jflUL8553b4qLHk2N2i4qBe9qQLNHgYM/ArtmAoN+4OEjMjOsaSKyYk556Wgft1Ltbw9kLZPZ000/cGwNEHsG1iAnJwfLli1Tm+wTWTIGTURWTLrl3HOTEescouqZyMwFNQLq9Qag0c7bZAUyMzPx6KOPqk32iWwqaNq2bRv69++P6tWrw87ODitXar/FlmTLli1o1aqVKgKsW7euWriRiCp/yZTO0doC3J0Bj0Bjx+9IFkE3/cA/P3JpFSIzU+q/oqmpqWjevDlmz55t0P3PnTuHfv36oUePHjh48CDGjBmDZ555BuvXc1gtUWWqn/y3dskUe0/s9+vLg20pavfQLq2SnQocWGTq1hBReQrB7733XrUZau7cuahVqxY+++wzdblRo0bYsWMHPv/8c/Tp06e0L09EBhaAd45Zqn7u9b+PS6ZY2tIqkm1a/TLw9zygw4uAA8fsEJmDSs/XR0REoGfPngWuk2BJrieiylE14yzqpexDHuwRwSVTLE/TRwH3ACDpMnBstalbQ0TGCpoiIyNRtWrVAtfJ5aSkJKSnpxf5GCkWlNvzb0RkOF0t0xGfO5HgHMxDZ2mcXIG2I7T7u78ydWuI6AazrAydMmUKfHx89FtYWJipm0RkMTxy4tEiYYPa3xnwqKmbQ2XV9hnAwRm4vFc72SURmVyld5QHBwcjKiqqwHVy2dvbG25ubkU+Zvz48Rg3bpz+smSaGDgRGVbL1D52FZw0Wbjk1ggX3DmZpcXyDAKaPqKd7HL3bCDMMkcdOzs745tvvtHvE1mySg+aOnbsiLVr1xa4buPGjer64sjUBFyjiKj0HPKy0D52hX6aAVVUTJZLCsIlaDq6Gki4BPhaXtbdyckJTz/9tKmbQWSa7rmUlBQ1dYBsuikFZP/ixYv6LNGQIUP093/++edx9uxZvP766zh+/Di++uortXjj2LFjK+YdEJFes8RN8M6JRaJjAA773sUjY+mCmwLhXQFNLrDn/0zdGiKbV+qgad++fWjZsqXahHSjyf7EiRPV5WvXrukDKCHTDfz+++8quyTzO8nUAwsWLOB0A0QVTaNB5+hlajci4GHk2XGYulXo+KL25/7vgMwUWBpZOkU+A2TjMipk6Ur9V7V79+7QaDTF3l7UbN/ymH/++af0rSMig4Wn/ouQjJPIsnPBXv8BPHLWol4fwL8OEHcGOLgYaP8sLImMhr7vvvv0PRWOjgzmyXKZ5eg5Iiq9Ljcms/zH7x6kOfrwEFpQUX9Jhf2wt7+5tMrfc4C8PKO1jYgKYshPZCFK+mD1y7yCRkk7bhaAk3VpPhjYPBmIOwucXAc05LI4RKbATBORFegc+wvsocEJr/aIdg03dXOoorl4Aq2Gavc52SWRyTDTRGThXHJT0Trud7W/M2CQqZtDFZhJnPJQvnm22j8HRMwGzm8Hrh0CqjXjsSYyMmaaiCxcm7g1cM1LQ5RLOE55tjV1c6iy+IQCje/X7jPbRGQSDJqILJidJhedYn5R+5zM0oamHzj8C5AcaerWENkcds8RWbDGSTvgn30NqQ4+atQcWbnQNkBoO+DyHmDvAuCuCTB3snTKl19+qd8nsmQMmogseNRc55if1c89VQYgx97FiK0ik+k4Cli2B9i3EOj6P8Cp6DU8zWkZlRdfvJEhI7Jw7J4jslDV046jVuq/yIUDIqo8bOrmkLE07A/4hAFpscBh7QzwRGQcDJqILFSXGO0H5iHfu5HsFGDq5pCxODgC7W7MCh7xlVo+x5zl5uZiy5YtapN9IkvGoInIAnllx6BZwp9qf0fAo6ZuDhlbqyGAkwcQfQw4u8Wsj39GRgZ69OihNtknsmQMmogsUIfY5XBALs65N8NV94ambg4Ze4kVN1+g5RPafU4/QGQ0DJqILIxjXibax65S+zsDOZmlzWr/vEw6AZzaAMScMnVriGwCgyYiC9Myfj08chMR51QNR727mLo5ZCpV6gD1b0wzsXsOzwORETBoIrIkGo1+moFdAQOhsXMwdYvIlDq8oP35709AWhzPBVElY9BEZEHqpexF1czzyLR3wz7/fqZuDplarTuBqk2A7DTgwCJTt4bI6jFoIrIguizTPv/7kOngaermkKnZ2d3MNu2ZD+Rmm7pFRFaNM4ITWYjAjAtokLwbebDDrioDTd0cMhdNBgJ/vgskXQGOrgKaDjS7GcE//fRT/T6RJWPQRGRGdMPLpzzUtNgs0zHvLohzCTF628h0ipp2QP874uQKtBkBbP1YWxBuZkGTrDf32muvmboZRBWC3XNEFsA9JxEt49ep/Z2czJIKazsCcHAGruwDLu3h8SGqJMw0EVlAZqFd3Co4azJxxbU+znm0MFm7yEx5BgFNHwEO/qid7DKsHcyFLJ1y4MABtd+qVSs4OHDEJ1kuZpqIzJxDXjY6xvyq9ncGPqot/iUqTFcQfnQ1kHDJbI6PLJ3Srl07tXEZFbJ0DJqIzFzTxM3wzolFkmMVHPK529TNIXMV3BQI7wpocoE980zdGiKrVKagafbs2QgPD4erqyvat2+PPXtK7kOfMWMGGjRoADc3N4SFhWHs2LH8xkFkCI0GXWKWqt2IKg8h156jj6gEHV/U/ty/CMhM4aEiMnVN09KlSzFu3DjMnTtXBUwSEPXp0wcnTpxAUFDQLfdfvHgx3nzzTSxcuBCdOnXCyZMn8fTTT8POzg7Tp0+vqPdBZLGKXJD1hvDUfxGSfhLZds7YU+V+o7aLLFC9PoB/bSDurHaW8HYjTd0iItvONEmgM3LkSAwbNgyNGzdWwZO7u7sKioqya9cudO7cGY8//rjKTvXu3RuDBw++bXaKiKDPMh3wuwdpjr48JFQg2NZtN/+i2wPtb9Q2yfQDeXk8YkSmCpqysrKwf/9+9OzZ8+YT2NuryxEREUU+RrJL8hhdkHT27FmsXbsWffv2LfZ1MjMzkZSUVGAjsjX+mVfQKGmH2uc0A2SwFo8DLj5A3Bng1AYeOCJTBU0xMTFq+GjVqlULXC+XIyMji3yMZJjef/99dOnSRc0GW6dOHXTv3h1vvfVWsa8zZcoU+Pj46DepgyKyNZ1if4E9NDjh1R7RruGmbg5ZChdPoPUQ7f7u2aZuDZFVqfTRc1u2bMFHH32Er776Ss3VsXz5cvz++++YPHlysY8ZP348EhMT9dulS+YzfJbIGFxyU9Ambo3aZ5aJSq3dc4CdA3BuGxD5n0kPoHxZnjRpktq4jArZVCF4QECAmpgsKiqqwPVyOTg4uMjHvPPOO3jqqafwzDPPqMtNmzZFamoqnn32Wbz99tuqe68wFxcXtRHZqrZxa+CSl44ol3Cc8jSfiQrJQviGAY36A0dXamubHpht0mVU3n33XZO9PpHJMk3yy9+6dWts2rRJf11eXp663LFjxyIfk5aWdktgpJsRVqPRlK3VRFbMXpODTjG/3MwycTJLKs/0A4d/BlKieQyJTNE9J9MNzJ8/H4sWLcKxY8fwwgsvqMyRjKYTQ4YMUd1rOv3798ecOXOwZMkSnDt3Dhs3blTZJ7me0+kT3apx4nb4ZUcixcEX//j14SGisgltC4S0BnKzgH1fm+woyhfrI0eOqE32iWxqnqZBgwYhOjoaEydOVMXfLVq0wLp16/TF4RcvXiyQWZowYYKak0l+XrlyBYGBgSpg+vDDDyv2nRBZic4xP6ufMi9Tjj27qen2dNMOTHmo6c0rJUPZYRTw6whg7wKgy1jA0fi/T+np6WjSpInaT0lJgYeHh9HbQFRR7DQW0EcmUw7IKDopCvf29jZ1c4gqVP55dkLTjuLF088ix84Rnzb8BclOATzaZLACQZPIzQZmNgeSrgD3fwW0fMLoR1N6Ijw9PdU+gyay9DiDa88RmcgtExPKZJbR2sksD/n2ZMBE5efgdHNWcCkIN//vyERmjUETkZnwyYpCk8Qtan9HwCBTN4esZZbwVkMBJ3cg6jBwXjtZKhGVDYMmIjPRMfZXOCAXZzxa4ppbPVM3h6yFuz/QfLB2f/dXpm4NkUVj0ERkBpxz09AubrXa38ksE1W0DjfWozvxBxB7hseXqIwYNBGZyWSWbrkpiHYOw3HvTqZuDlmbgHpAPZm+QgNEcGkVIqNNOUBEFT+ZZecYbQH49sDHoLHjdxmqBJ1eBk6tBw7+CPR4C/AwzshMWTrl1Vdf1e8TWTIGTUQm1jThL/hlR92YzPIeUzeHrFV4F6B6S+DqP8Ce+UCPm5MQVyZZSWLq1KlGeS2iysavtESmpNGga/RPajciYCAns6TKI5NddnpFu79nHpCVxqNNVEoMmohMqE7KfoRknESWnSt2V3mQ54IqV6MBgG9NID1O201nBLJ0yvnz59XGZVTI0rF7jsiICk9meWf0YvVzn38/pDn68FxQ5XJw1C7k+8fr2oLwNsMBe+0C6pW5jEqtWrXUPmcEJ0vHTBORiQSnn0b9lD3Igz12BDzK80CV4pbJLls+Cbj5AfHngGO/8agTlQKDJiIT6RqzRP38z6cb4l1CeB7IOJw9gLbPaPd3fcGlVYhKgd1zRCbgnXUdzeM3qv3tgY/zHFCly59tmtL7WWDnF8CV/cDFCKAm5wYjMgQzTUQm0Cn2F7VkylmPFrjs3ojngIzLMwhocWNpFQmeiMggzDQRGZlLbirax65S+9uYZSIjDTq4RceXgf2LgJN/ANEngMAGPBdEt8FME5GRtYtbBde8VES5hOOkVwcefzKNgLpAw37a/V2zeBaIDMBME5ERv/U75GWjc8wytc8lU8jkZLLL42uAQ0uBuyYAXsEV/hKOjo4YNWqUfp/IkvE3mMiIWiSsh092NJIcq+Cgb28eezKtGu2BsPbApb+B3XOAXu9V+Eu4uLhg9mwuEkzWgd1zREaaI8dOk4tu13/UZ5ly7Z157Mn0uozV/tz7NZAeb+rWEJk1Bk1ERnJH4lYEZl1CmoMX9vjfz+NO5qFeHyDoDiArGdizoMKfXqPRIDo6Wm2yT2TJGDQRGYNGgx7Xv1e7u6oMRJaDO487mQd7e6DrOO3+33OArNQKffq0tDQEBQWpTfaJLBmDJiIjqJ/yN6pnnEKmvRsiAgbymJN5afwA4FcLSIsFDnxn6tYQWVfQJEV94eHhcHV1Rfv27bFnz54S75+QkIAXX3wR1apVU0WB9evXx9q1a8vaZiKL0/1GlmmP/wAuzEvmuZBvlzE3J7vMyTJ1i4isI2haunQpxo0bh0mTJuHAgQNo3rw5+vTpg+vXrxd5/6ysLPTq1Qvnz5/HL7/8ghMnTmD+/PkICeFaW2Qbaqb+i1qp/yLHzhE7Ah8zdXOIitZ8MOBVDUi+ChzSrotIROUMmqZPn46RI0di2LBhaNy4MebOnQt3d3csXLiwyPvL9XFxcVi5ciU6d+6sMlTdunVTwRaRLeh+/Qf184DfvUhyCjR1c4iK5ugCdHxJu79jBpCXyyNFVJ6gSbJG+/fvR8+ePW8+gb29uhwREVHkY1avXo2OHTuq7rmqVauiSZMm+Oijj5Cby/+QZN3TDMhWLf0UGiZHIA/22Bb4hKmbRlSy1k8Dbn5A3BngqHapHyIqY9AUExOjgh0JfvKTy5GRkUU+5uzZs6pbTh4ndUzvvPMOPvvsM3zwwQfFvk5mZiaSkpIKbESWqNuNLNNhnx6IdQk1dXOISubiCbR/Qbu/fboa9UlERpwRPC8vTw01nTdvHhwcHNC6dWtcuXIFU6dOVXVRRZkyZQree6/iZ6YlMqbAjAtomrhZ7W8JepIHn8x6Ud8pDzXV7rQbCez6Aog6DJxcDzS4p1yvIUunDB06VL9PZMlK9RscEBCgAp+oqKgC18vl4OCi1yySEXNOTk7qcTqNGjVSmSnp7nN2vnVW5PHjx6ticx3JNIWFhZWmqUQm1+P6IthDgyPeXRHpVs/UzSEqMli6hbs/0HYEsHMmsPVjoH4fwM6uzEdPRkx/++23PPpke91zEuBIpmjTpk0FMklyWeqWiiLF36dPn1b30zl58qQKpooKmHT/yby9vQtsRJYkIOMimif8qfY3Bz1t6uYQlX4hXyd34Oo/wKmNPHpEZR09JxkgmTJg0aJFOHbsGF544QWkpqaq0XRiyJAhKlOkI7fL6LnRo0erYOn3339XheBSGE5kre5SWaY8HPXqjKvuDUzdHKLS8QgA2j6j3ZdsUzlqm2TpFPmMkI3LqJClK3UH86BBg9QaQhMnTlRdbC1atMC6dev0xeEXL15UI+p0pFtt/fr1GDt2LJo1a6bmZ5IA6o033qjYd0JkJgIyJcuk/Xa+qar2ywSRRWab9i4AruwHTv8J1OtVpqeRpVM8PT3VfkpKCjw8PCq4oUTGY6exgNBfapp8fHyQmJjIrjoyewdmDkKr+HU45tUJ39X61NTNITKIvhA8vw0TgF2zgJA2wDN/lqm2STJMDJrIWuIMrj1HVJFiz6B5vC7LNJzHliw/2+ToBlzZB5y5WctKZKsYNBFV4IikfT9MgANycdyrI664N+SxJcvmGaQdSSe2lK+2icgaMGgiqiD+mVfQMn692mctE1ka3Sz2xWabLu8FzmjnHSOyVQyaiCpwXibJMp3w6oDL7o15XMnilwFSvKoCbW50NTPbRDaOQRNRRYg5dTPLxHmZyNp0Hg04ugKX9wCnWdtEtotBE1FF+OsjlWWSeZkueTThMSXrItkm3bxNm96TWY0NfqisBjFw4EC15V8ZgsgScSEgonKQLoxq6SfxyqnlyIMdNgaP5PEk69RlHLB/ERB5CDi2CrjjQYMe5urqimXLllV684iMgZkmonLqHTlf/Tzk2xORbnV5PMlqFKht8qgCdHpJu7/5QyA3x6RtIzIFBk1E5VAz9V80TI5ALhzwZ9UbQ7OJrFXHFwH3KkDsKeDfn0zdGiKjY9BEVFYaDfpEzlO7+/z7IdYllMeSrJuLl7abTjeSLjvDoBnB7ezs1Cb7RJaMQRNRGUiXxTffLUCt1H+RbeeMzRwxR7ZCCsK9Q4Cky8C+haZuDZFRMWgiKgM7TR56X9NmmSKqPIQk5yAeR7IJ4387heVeT2gvbP8MyEw2dZOIjIZBE1EZ3JG4FSEZJ5Fp74atQU/yGJJN2e/fF/CvA6TFALvnmLo5REbDoImotHKy0Cfy/9Tu9oDHkOboy2NINiXPzhHo8Zb2wq5ZQGqMqZtEZBQMmohKa99CBGRdRrKjH3YEPsbjR7bpjoeA4GZAZpK2KJzIBnByS6LSSI8Htmo/IGSKgUwHDx4/snpFLuRrbw/0+RBY1F9bEN7uWSCwvimaR2Q0zDQRlYYUvqbHI8olHPv87+OxI9tW606g/r2AJhfYOLHIu8jSKX379lUbl1EhS2en0Wg0MHNJSUnw8fFBYmIivL29Td0cslXx54Ev2wK5WfgmfCpOenc0dYuITG7Kna7AVx2AvBxgyGqgdjdTN4mo0uIMZpqIDPXnuypgQu0eOOnVgceNSATUA9oM1x6LDW8Debk8LmS1GDQRGeLSHuDICrUo70yHoYCdHY8bkU63NwEXHyDyMPDvEh4XsloMmohuR3qw12uHV+/368tFeYkKk8V87/yfdn/zZCDr5nIpsnSKh4eH2riMClk6Bk1Et3NkBXB5L7LsXLExeCSPF1FR2j0H+NYEkq8Bu74scFNaWpraiGwyaJo9ezbCw8Ph6uqK9u3bY8+ePQY9bsmSJWrRxgceeKAsL0tkfFlp+lFB24IeR7JTAM8CUVGcXIGe72r3d84AEi/zOJHVKXXQtHTpUowbNw6TJk3CgQMH0Lx5c/Tp0wfXr18v8XHnz5/Hq6++iq5du5anvUTGteNzIPES4BOGbYGP8+gTFTGHk27DHQ8CNToC2WnAhgk8VmR1Sh00TZ8+HSNHjsSwYcPQuHFjzJ07F+7u7li4sPjVrnNzc/HEE0/gvffeQ+3atcvbZqJKJx8AU39ah+ztM9TlH3yfR7a9K488UUlkgETfqYCdvbZb++xWHi+y3aApKysL+/fvR8+ePW8+gb29uhwREVHs495//30EBQVhxIgR5WstkRH1u/oFnDRZOOXZFke87+SxJ7oNlXHaBaDtM9or/ngdyM3mcSPbXEYlJiZGZY2qVq1a4Hq5fPz48SIfs2PHDnz99dc4ePCgwa+TmZmptvyTThEZU4OkXWicvBO5cMBv1cdwigGi0pDFfP/7FYg+Duz9mseOrEaljp5LTk7GU089hfnz5yMgwPAC2ilTpqiZOXVbWFhYZTaTqKCcTNx39Qu1uyNwEKJda/IIEZXC+D8u41d/7UhT+x3T0a1LR3Tr1k31TBDZTKZJAh9ZOygqKqrA9XI5ODj4lvufOXNGFYD3799ff11eXp72hR0dceLECdSpU+eWx40fP14Vm+fPNDFwIqPZ8TkCsi4jybEKNgc9zQNPVAYyp1m72FUISz+GLWMaAw8v4HEki1eqsN/Z2RmtW7fGpk2bCgRBcrljx1vX4WrYsCEOHz6suuZ024ABA9CjRw+1X1wg5OLiotZ+yb8RGUX0Se2ivADWVH8FWQ7uPPBEZaCxs8eqkP9pi8IPLwNO/8njSLaVaRKSARo6dCjatGmDdu3aYcaMGWqWVxlNJ4YMGYKQkBDVxSbzODVp0qTA4319fdXPwtcTmZxGg7PfjkTt3Cyc8OqAwz53mbpFRBbtintD7aSXf88B1owDRu0GnPlFhGwoaBo0aBCio6MxceJEREZGokWLFli3bp2+OPzixYvstybL9M8PqJ16UM38vVJ9Q+b6ckTlldp+DMIfmQZo/sP52u/DY8DHPKhksew0GllYy7xJTZMUhCcmJrKrjipHSjTwZRsgIwFrq43Cdk5kSVQhJvSpDU9PT+1/s7d84fHyNiC4KY8uWWScwaEMREIW5M1IwFXXetgZ8CiPCVFl0OQCq18B8nJ5fMkiMWgiOvEHcPhnSbxiRejryLMrda81ERnCxQu4egDYPYfHiywSgyaybenxwG9jtPudXsJl90ambhGRVZm46r+bF+7WLn6NzZOBmFMmaxNRWTFoItv2x5tASiRQpR7Q421Tt4bIurV4AqhzF5CTAax8gd10ZHEYNJFtd8sdWqKdR+aBOYCTm6lbRGTdZETqgFmAizdweS8Q8aWpW0RUKgyayDalxQG/jdbud3wJCGtr6hYRWSU7O3sE17lDbZNWH8X4TfH4JXCU9sbNHwLRJ0zdRCKDseKVbI/MsrH2NSAlCgioz245okrk5OKKpz/9qcB1+/36oWniFjRI3g2seB4YsQFwcOJ5ILPHTBPZnkNLgf9+AewcbnTLuWL88sNqIyIjsLPD8tA3AFcf7Wi6rZ/wsJNFYNBEtiXuHPD7q9r97uOB0DambhGRTUpyCgTum6G9IOs9nt9p6iYR3Ra758h25GYDy0cCWclAjY54K6YXNMwuEVWq7Mx0zB/9oNofOXMFnFxuDrgYf7IeHvbrizbxa4HlzwIv7ADc/HhGyGwx00S2Y9tUNWIn3d4Tn7j/DxrpniOiSiUrdSVFX1VbUat2/VZ9DGKcQ4Gky8CasdqaQyIzxaCJbIOk/iVoArAy9FUkOAebukVEBCDLwR1La0wC7B2BIyvUwtlE5opBE1m/5Cjgl2GAJg9o/jgO+fY0dYuIKB81E79uctm1rwKR+WYRJzIjDJrIuuXmAL+O0E4vENgI6DfN1C0ioqJ0HgPU7amdLfznp4CMRB4nMjsMmsi6/fUhcH474OwJDPoecPYwdYuIqCj29sBD8wGfMCDuLP776imM//UQpwIhs8KgiazXiXXAjunafVm6IaCeqVtERMVQc6WtuwI8sgiwd0KTpK3oErOEx4vMCoMmsk4xp4EVz2r32z8PNHnI1C0iskl2dnYICK2jNtm/nfF7nLEq+GW1f8+1uaiVcsAIrSQyDOdpIuuTngD89Ji2JiK0HdBrsqlbRGSzZF6mZ2auKNVjdld5EDXS/kPLhA144sI7QNydgH+tSmsjkaEYNJH1FX7/MhyIPQV4hwKDfgAcnVkXQWSBy6wEZl5EaPpx7ZegERsBV29Tt4xsHLvnyLpsnAic2QQ4uQODFwNeVU3dIiIqgxx7F3wX/jESHQOA6OPAr88Aebk8lmRSDJrIehz4Dtg9W7svC/FWa27qFhHZPFlGZcHoB9Um+6WR7BSAH8KnINvOGTi1Hvhzks0fTzItds+RdTi5AfhtzM2FeO94QD8ih4hMR5ZOibl8Rr9flokvfwl7C4MvvgvsmgX41ADa3xjkQWRkDJrI8l3eDywbCmhygeaDMT62L8BgichqyCz+AZmX0StqAfDH64BnkP6LEZHZd8/Nnj0b4eHhcHV1Rfv27bFnz55i7zt//nx07doVfn5+auvZs2eJ9ycqldgzwOJHgOw0oM7d2vmYDBjWTESWZXPQUOz2l0BJg5xfngHObb85v9ONjcjsgqalS5di3LhxmDRpEg4cOIDmzZujT58+uH79epH337JlCwYPHoy//voLERERCAsLQ+/evXHlypWKaD/Z+ppyPzwEpMUC1VoAj34HODiZulVEVBns7LA6ZCyOeHeFoyYbWPI416gj8w+apk+fjpEjR2LYsGFo3Lgx5s6dC3d3dyxcuLDI+//4448YNWoUWrRogYYNG2LBggXIy8vDpk2bKqL9ZKtSY4DvBgDx5wG/cOCJZYCLp6lbRUSVSGPngCU13sU592ZAZhKSF/RHQMZFHnMyz6ApKysL+/fvV11s+iewt1eXJYtkiLS0NGRnZ8Pf37/Y+2RmZiIpKanARnTzlygO+O4B7TBkr+rAUyu0NQ5EZDNTEVx1rQuvnDg8c/YVVMm8bOpmkY0oVdAUExOD3NxcVK1acO4buRwZGWnQc7zxxhuoXr16gcCrsClTpsDHx0e/SZcekSKzfEuXXNRhwCMIGLoa8K/Ng0NkpmTpFO/A6mozZBkVQ2Q4emNh7c8R6VILPjkxKnDyzbrG2iayrnmaPv74YyxZsgQrVqxQReTFGT9+PBITE/XbpUuXjNlMMuflUX54GLj6D+BeRRsw5VuEl8WgROa5jMqouevUJvsVJdXRD1/XnonrLjXhm30dI89I4GTYl3ciowRNAQEBcHBwQFRUVIHr5XJwcHCJj502bZoKmjZs2IBmzZqVeF8XFxd4e3sX2MjGSQ3TovuAy3sBV1/gqZVAUCNTt4qITCjFyR8Las9EjHMo/LOv4bkzo1jjROYTNDk7O6N169YFirh1Rd0dO3Ys9nGffvopJk+ejHXr1qFNmzblazHZnsQrwDf3ApHSJRcIPL0GqFZy4E1EtkFmDZ9fZxaiXMJVxkkCpy9+XM7MM5nH5JYy3cDQoUNV8NOuXTvMmDEDqampajSdGDJkCEJCQlRdkvjkk08wceJELF68WM3tpKt98vT0VBtRieLOAt/dDyRcBLxDgCGrbumSIyLzlZ2ZgR/f0X4+PDH5Gzi5FF+aUVZJToEqcBp27n8IST+JkWdexqJan+KCB79ckYlrmgYNGqS62iQQkmkEDh48qDJIuuLwixcv4tq1a/r7z5kzR426GzhwIKpVq6bf5DmISnRpL7CglzZgkmLv4esKBExEZP40mjxEnjmiNtmvLFLjNL/2Fzjn0RxueSkYfnYsGidurbTXI9tkpynLYkBGJlMOyCg6KQpnfZONOLoKWP4skJMBBDfTzsPkdWvdHDNNROYtKyMN05/ooPbH/bgbzq7ulfp6TnkZePzCO2iYHIE82OGPaqOwI+AxTHmYWScqf5zBtefIvEgML4tybpyolktAvT7AwIW3TFzJYImIipJt74rvw6fgvqsz0TF2Bfpdm63WrXv717HIsyv4kTfloaY8iFQqDJrIfGSnA2vGAf8u1l5uOxK452PAQftrykCJiAwhwdHq6uMQ6xyGvtdmoX3cKlTJuqxmE5duPCKLmKeJqFiyHMrXvbQBk5090GcK0HeqPmAiIioVOzvsDHwUP9T8CFl2rqibsh8vnRqB0LSjPJBUZgyayPRObQT+r5t2SgH3AO0IuY6j1B89IqLyOObTFbPrzUO0c9iNKQleRLvYldpSAKJSYtBEppOTCax/G/hxIJCRAIS0AZ7bBtS6k2eFyIq4efupzVSuu9bG7HoL8J/3nXDUZOPBK9Mw+OIkvPfzTnb7U6lw9ByZxvVjwK/PAFH/3axf6vMh4Ohyy11Zy0REFUKjwZ3Ri9E7ch4ckIsEpyAsC5uAs56tbrkri8RtS5KBo+eYaSLjyssFIr4C5nXXBkzSHTd4KdBvWpEBExFRhbGzw7agJzC37hy19Ip01404Oxr3Xp0Nx7xMHmi6LVbZkvFEHQFWvwxc2a+9XLcX8MBXgGdQkXdnhomIKsNl98aYVW8h+l2bhXZxv+HOmJ/QOGk7VoS+hrOerYv9W8TsEzHTRMaZSmDzB8D/3akNmFy8gfs+105YWUzARERWtIzKxOFqk31zkeXgjhWhb+C7mh8j0TEAAVmXMfLsaDx86SO45SSZunlkpphposojo1NkZu8N7wCJF7XXNein7Yrzrs4jT2QDZOmUS0f26ffNzTGfLjjr2QL3RM5Fh9iVaBO/Fo2TdmDVvBHYU+X+AhNi5s9+M+tkm5hpospx7RDw7X3AsqHagEkW231kEfDYjwyYiMisZDp4YlXIq5hb5ytEutSCe24S7r/6OV4+OQx1k/cW+RgJoFhCYHuYaaKKFXMK2DIF+G+5dhkUR1eg8xig8yuAs0eBu/IPDhGZkwsezTCr/jdoF7sKvaK+RnDmOYw4NxanPNtgY9VncMmjiambSCbGoIkqRtxZYNs04N+fAF0KvslAoOe7gG8YjzIRWQTpjtsd8DAO+vZCz+vfoEPMctRL2ae2414dsbHqCFx1b1jilz923VkvBk1UPlcOALu+0NYu6YKlBn2BHm8BwVwMk4gsU4ajN9ZUH42dAY/grqhv0TJ+PRomR6hNMk/bAx/HKc+2XLnAxjBootLLzQFOrgP2/B9wbtvN62UKge7jgdBbh+zmx245IrIU8c7V8WvYW9gS9BTujvoWzRL+1GeerrrWxc6AR3HY9y5k27ve8lgWjlsfzghOhku6Bhz4DjiwCEi6or3O3lHbDdfpZSDYsP5+Bk1EtiMrIw2zhndX+y8v3AJnV3dYMt+sa+gS8zPaxK2BS166ui7d3hP/+PXBXv/+iHSrW+Lj2XVn2TOCM2iikmUkAcd+Aw7/rM0q6brg3KsALZ/ULn9iQM0SAyUisiYyl1O7uFVoF7ca/lnX9NdfdG+MA3734j+f7kh19DMoaGJGynKCJnbP0a0yU4Azm4H/ftV2w+Xkm5CuRkegzQig8QAue0JENivd0Rtbg57CtsAnUCdlvwqgGiduR420o2obcOVznPFshcM+d+GITzekOfqoxzFAsmzMNJFWwkXg5HrgxB/A+e1AbtbNIxNQH2j6KNB0IOBf67ZHjFklIrJFntlxaBm/Ds0SNyE0/YT++lw44IJHU5zw6oATXh0R5Vq72AJydt+ZBrvnqGSpMcD5HTe27UD08YK3+9UCGvYDmj4CVGtu0AgRBktEVFhOViZWTB2n9h98bTocnW1jYW7/zCtomrgZzRI2o3rGqQK3JToF4qRne5z1bIVzHs2R6Fz1lsczeDIuds/RTXm5QPQJ7bpvVw8AF3cD148WPEJ29kBYB6DBPUD9e7TZJQMCJcFgiYiKk5eXizMHtuv3bUWcS4jqvpNNAqgGybtRP3k3aqccgE92NNrGr1GbiHcKxnmP5jjn0QyX3Rvhukstg/++MrgyLtY0WWPhtmSNJCi6fhyIPAxc/QfITr31vkF3AOFdgFpdgZqdAXd/U7SYiMjqA6gIl4cREfAwHPMyUSv1IOol70F46iFUTz8Jv+xI+CVEomXCenX/HDsnXD5TG1fcGqCdW31cd62F6y41keboW6rXZf2UmQRNs2fPxtSpUxEZGYnmzZtj1qxZaNeuXbH3X7ZsGd555x2cP38e9erVwyeffIK+ffuWp922LTsdiL8AxJ+/ucWe0gZJSZeLfoyzJ1CtBRDSCghtow2SPAJK/dLMKhERlV2OvQtOebVXm/rTnJuGGmlHEJ76L2qmHUZI+gm45aaomqj8dVEixcEX0a41cd0lHNEuNRHnXB0zfzyNeOdqePeRDiX+ndZdlz8zVfg6BlmVEDQtXboU48aNw9y5c9G+fXvMmDEDffr0wYkTJxAUFHTL/Xft2oXBgwdjypQpuO+++7B48WI88MADOHDgAJo04To+BUjqOj0BSIkCkq8ByZHan/kvJ1wCUiJLPkle1YGghkBQYyCoERDSWtvdZu9w2/PLoIiIyHiyHNxx2qut2hSNBn5ZVxGSflIFUNUyTiEo44LKRnnmJsAzNQG1Uv+95XnSjnsj3jlYdfX1dwpEsqM/UpyqqJ9qcwpAqqNviQEVVcLoOQmU2rZtiy+//FJdzsvLQ1hYGF5++WW8+eabt9x/0KBBSE1NxZo12r5b0aFDB7Ro0UIFXhVZoGUycggl+5OVCmSlANlpN/bzbdk3fkr3WXp80VtGonaRW0O4eAN+4Tc3/9raACmwIeBmWAq3qG8V/M9DRBU9ueX0J7RZkHE/7rb4yS1NxSkvHYEZFxGUeQFBmecQkHkJflnX4JcVCY9c+ey4vTzYIcPBE+kOXsVs3kh38ESWvVu+zVX/8/X+rQAnd7y15gw0UgdrRbVVlVIInpWVhf3792P8+PH66+zt7dGzZ09EREQU+Ri5XjJT+UlmauXKlTAbu74Erh/TDrNXWzaQm3lzPyfzxnVZt245WTfmMSpV7FkyN3/AqxrgFZxvqwZ4VgV8QrQj29z8KnREG4MlIiLzlW3vhqvuDdRWmHTxqbqoG0GUZ04svLLj4JUjW4yaCsEzJx4OyIV7brLayuSE9sdHEg/YuSDX3knVX+XaOSH2o5v7OfY3ftoV/Jln56CCrTzYo02tAO2KEnYO2l4Q9dO+6OvUT0egWjOgtnZ2eVMpVdAUExOD3NxcVK1acHikXD5+vNCQ9Ruk7qmo+8v1xcnMzFSbjkR+ukiwUhz6Hbiws2Key9ENcPYAnNwAJw/A2V27OXlqr3Px1AY8kg1y9dXuu+a/7AM4OKmnenf1ESBFOrIB6CeclfmTTuDdAXfoX1Ldj4jIDGVlaJcaEZlpqdDk3VhVgCqMfFomIwgXHYMAx+ZF3sdOkwe33CS45yTBNS9Z1U255qbc+Hnjcp5cToZTXgac89LhnJepMlxOeZlwzsuAiybfRMeQfe1l+xub9pPLMEm3qTIpUqshQEArVAZdfHG7zjezHD0n9U/vvffeLddLN6D5K2MEX0qfG+VViIgqzuyRd/NwUjnMvrFVnuTkZNVNVyFBU0BAABwcHBAVFVXgerkcHBxc5GPk+tLcX0j3X/4uPambiouLQ5UqVWBn4NxBZDwSoUtAe+nSJfOsOaMy47m1Xjy31ovntvQkwyQBU/Xq1Uu8X6mCJmdnZ7Ru3RqbNm1SI+B0AY1cfumll4p8TMeOHdXtY8aM0V+3ceNGdX1xXFxc1Jafr2/p5qcg45OAiUGTdeK5tV48t9aL57Z0Ssowlbl7TjJAQ4cORZs2bdTcTDLlgIyOGzZsmLp9yJAhCAkJUV1sYvTo0ejWrRs+++wz9OvXD0uWLMG+ffswb9680r40ERERkcmUOmiSKQSio6MxceJEVcwtUwesW7dOX+x98eJFNaJOp1OnTmpupgkTJuCtt95Sk1vKyDnO0URERESWpEyF4NIVV1x33JYtW2657pFHHlEbWSfpSp00adItXapk+XhurRfPrfXiuTWjyS2JiIiIbFHBKT2JiIiIqEgMmoiIiIgMwKCJiIiIyAAMmqhE27ZtQ//+/dWEXzKx6O3WDJSBAHK/wltJy+aQaci0ILL4tpeXF4KCgtTcaydO3FhcqgTLli1Dw4YN4erqiqZNm2Lt2rVGaS9V7rn99ttvb/l/K+eYzMucOXPQrFkz/RxMMufhH3/8UeJj+H+24jBoohLJHFzNmzfH7Nmlm7pe/kBfu3ZNv8kfbjIvW7duxYsvvojdu3erCWezs7PRu3dvdc6Ls2vXLgwePBgjRozAP//8oz6MZfvvv/+M2naq+HMr5EM4///bCxcu8FCbmdDQUHz88cfYv3+/mvPwrrvuwv33348jR4peg5T/ZysWR8+R4b8sdnZYsWKFfjb44jJNPXr0QHx8PGdxtzAy/5oEt/KBe+eddxY7T5t88K5Zs0Z/XYcOHdR8bXPnzjVia6miz61kmmTlhoSEBB5cC+Pv74+pU6eqLzOF8f9sxWKmiSqFfIhWq1YNvXr1ws6dO3mULUBiYqL+D3BxIiIi0LNnzwLX9enTR11Pln1uRUpKCmrWrKnWkiwpe0HmITc3V62yIV9kiluajP9nKxaDJqpQEihJxuHXX39Vm/zx7d69Ow4cOMAjbcZkDUnJMnTu3LnE2fqlNk03+7+OXGbNmuWf2wYNGmDhwoVYtWoVfvjhB/U4WdHh8uXLRm0v3d7hw4fh6empJrF8/vnnVQ9A48aNi7wv/8+awYzgRCX94ZVNR/7onjlzBp9//jm+//57HjgzJfUvUpe0Y8cOUzeFTHRuJVORP1sh/3cbNWqE//u//8PkyZN5XsyI/I09ePCgyiD+8ssvaj1Y6XotLnCiisNME1U6Wdj59OnTPNJmSpZEkhqlv/76SxWZliQ4OBhRUVEFrpPLcj1Z9rktzMnJCS1btuT/XTPk7OyMunXronXr1mqkpAzWmTlzZpH35f/ZisWgiSqdfCOSbjsyL7KCknyoSmp/8+bNqFWr1m0fI5mITZs2FbhORmcVV09BlnNui6qXkW4g/t81f9KVmpmZWeRt/D9bsdg9R7ctDM2fJTp37pwKgqSgtEaNGhg/fjyuXLmC7777Tt0+Y8YM9Qf6jjvuQEZGBhYsWKD+aG/YsIFH2gy7bRYvXqxqWGQ+H11dko+PD9zc3NT+kCFDEBISor7NitGjR6Nbt2747LPP0K9fP1WEKsOe582bZ9L3QuU/t++//74aCSkZDBlBJ6OxZMqBZ555hofXjMjf3HvvvVf9/U1OTlbnWUYtr1+/Xt3O/7OVTBbsJSrOX3/9JQs637INHTpU3S4/u3Xrpr//J598oqlTp47G1dVV4+/vr+nevbtm8+bNPMBmqKjzKts333yjv4+cW9251vn555819evX1zg7O2vuuOMOze+//26C1lNFn9sxY8ZoatSooc5r1apVNX379tUcOHCAB9rMDB8+XFOzZk11ngIDAzV33323ZsOGDfrb+X+2cnGeJiIiIiIDsKaJiIiIyAAMmoiIiIgMwKCJiIiIyAAMmoiIiIgMwKCJiIiIyAAMmoiIiIgMwKCJiIiIyAAMmoiIiIgMwKCJiIiIyAAMmoiIiIgMwKCJyEJ1794dY8aMMXUzzEJsbCyCgoJw/vx5dZnHpniPPfaYWnCZiEqPQRORhVq+fDkmT54Mc2dIAJObm4tOnTrhoYceKnB9YmIiwsLC8Pbbb5f4+A8//BD3338/wsPDUZkiIyPx8ssvo3bt2nBxcVFt69+/PzZt2lTgfpcuXcLw4cNRvXp1ODs7o2bNmhg9erQK7kxtwoQJ6njJsSWi0mHQRGSh/P394eXlBWvg4OCAb7/9FuvWrcOPP/6ov14CFHmfkyZNKvaxaWlp+PrrrzFixIhKbaNksVq3bo3Nmzdj6tSpOHz4sGpvjx498OKLL+rvd/bsWbRp0wanTp3CTz/9hNOnT2Pu3LkqsOrYsSPi4uJgSk2aNEGdOnXwww8/mLQdRBZJQ0QWqVu3bprRo0fr91966SV12dfXVxMUFKSZN2+eJiUlRfP0009rPD09NXXq1NGsXbv2lud48cUX1ebt7a2pUqWKZsKECZq8vDz9ff744w9N586dNT4+Php/f39Nv379NKdPn9bfnpubq/nkk0/U8zs7O2vCwsI0H3zwgbpt6NChGvkzk387d+5cse9p5syZGj8/P83Vq1c1K1eu1Dg5OWkOHjxY4nFYtmyZJjAwsNhjI9asWaPe3w8//KAuJyUlaR5//HGNu7u7Jjg4WDN9+vRbHlPYvffeqwkJCVHHtLD4+Hj9/j333KMJDQ3VpKWlFbjPtWvX1Os9//zzmrKSY/3RRx9pwsPDNa6urppmzZqp968jbXvqqac0Hh4e6n1NmzatyPf13nvvabp06VLmdhDZKmaaiKzEokWLEBAQgD179qgMzQsvvIBHHnlEdXsdOHAAvXv3xlNPPaUyM4Uf5+joqB43c+ZMTJ8+HQsWLNDfnpqainHjxmHfvn0qW2Jvb48HH3wQeXl56vbx48fj448/xjvvvIOjR49i8eLFqFq1qrpNnk+yKyNHjsS1a9fUJl1axZF2N2/eXLXz2WefxcSJE9Xlkmzfvl1lgIoj7Rk8eLDKYD3xxBPqOnk/O3fuxOrVq7Fx40b1HHKMiiPZIckqSUbJw8Pjltt9fX3191u/fj1GjRoFNze3AvcJDg5Wr7906VL5soqymDJlCr777juVuTpy5AjGjh2LJ598Elu3blW3v/baa2p/1apV2LBhA7Zs2VLk+2rXrp0635mZmWVqB5HNMnXURkQVk2nKnznIyclR2QbJOuTPdMh/+YiIiALP0ahRowKZpTfeeENdV5zo6Gj1PIcPH1YZGxcXF838+fMNaqchjh07pp6/adOmmuzs7Nve//7779cMHz68yNf88ssvVYZsy5Yt+tukzZLByp+hSUhIUFmg4tr5999/qzYtX768xLbs3r1b3W/FihVF3i4ZLbk9KipKc/HiRf3xl/f6888/l/jcGRkZqo27du0qcP2IESM0gwcP1iQnJ6tMX/7niY2N1bi5ud3yvv7991/VjvPnz5f4mkRUkKOpgzYiqhjNmjUrUCNUpUoVNG3aVH+dLvtz/fr1Ao/r0KED7Ozs9JclMySjq6Q4W55HanMk4/P3338jJiZGn2G6ePGiylpJtuLuu++usNO4cOFCuLu749y5c7h8+fJti7vT09Ph6up6y/W//PKLeq+SUWrbtm2BmqPs7GyVbdHx8fFBgwYNin2N0maGDLm/ZPdmzJiBFi1aqAJzyZb17dtXZbIkK/bcc8/p7/vHH3+o2i453r169SrwPFlZWWjZsiXOnDmj9tu3b6+/TR5T1PvSZcEKZx2JqGQMmoishJOTU4HLEgjlv04XGOmCHkPJ6DAZ/TV//nw1GkweL8XE8gFduAuqvHbt2oXPP/9cdS198MEHqrj7zz//LBDUFSZdkvHx8bdcL4GEdE1JECaF2SU9x+3Uq1dPPf748eMl3q9u3brqfseOHVNdmIXJ9X5+fggMDFT3q1atmr7rTt6HdO9J0DRgwIACwU9ISAgOHTqk9n///Xd1OT8ZyVeaAnPdfaUdRGQ41jQR2TjJIOW3e/duFSRIlkmGyJ84cUINU5dsUqNGjQoEKHI/CZwKD7nPT4bcS9bqdiTr8fTTT6taLBmRJiPipO5G6ndKIsGR1FIVJiPE/vrrL1XfI7VSOjJdgASTe/fu1V8nw+9PnjxZ7GtIxqZPnz6YPXu2qvEqLCEhQf2U7J5kgr766iuVActPskmSQRo0aNAtAdz+/fvVMdLVe8moSAnAdJsc48aNG6vgSDJ8+W+TTR4n71feV/7zKeeqqPf133//ITQ0VAVqRGQ4Bk1ENk4+hKUwWoIjGSI/a9YsNaeQkKyIBALz5s1TQ+dluL3cV0e6xd544w28/vrrqkBZuogk6JKAR0e61+SDXIbs5+/eK0wKyqVbS4rKdY+bNm2aem7dpJVFkWBGiqKLyjbVr19fBU6//vqrfq4oCUiGDh2qiqblNnmsZLSkwL2kbJQETBLYSLeePJ90W0rm6IsvvlBdmjpffvml6rKUdm3btk3N2SRF5BJMSYZI5kgqnPUZMmSIOsYlkXa/+uqrqvhbivflWEsmTc6XXPb09FTvQ96XnCcJjCQIlfdVmBS+y8AAIiqlQjVORGShheCFi31r1qyp+fzzzwtcV7hIWR43atQoNQxehuTLcP+33nqrQGH4xo0bVbGyFHzLEHcpqs7/PDIMXqYYkNeTAusaNWqoYfE6J06c0HTo0EEVJBc35YA8p4ODg2b79u233Na7d2/NXXfdVaBNhbVr104zd+7cIo+NOHr0qJqGYdy4ccVOOSDP8eabb2pKIlMhyPQM8l6l6FqmIBgwYIDmr7/+KnA/KbCW6RaqVq2qjolMw/Dyyy9rYmJibinu7tq1q+a7777TGEKOwYwZMzQNGjRQzytTLfTp00ezdetWdbsUgz/55JPqfclrf/rpp7cci/T0dFUcn39AABEZxk7+KW2gRUTWQWbrlkJkKUi2ZFLnIxkWya4UlVm5HelykyyQFMBX9iSZOvKn9/HHH1eF2u+++67RzvGcOXOwYsUKVTdGRKXDQnAisnj9+vVT3WVXrlwpcR4onX/++UcVdUtXm9Qzvf/+++p6WYrFWGRUn8zZJKMeV65cqa77/vvvC4x4rAxS9yRdekRUegyaiMgqlHbxYqmXkjouKVSX4f5S52PMwuguXbqUeiRjRXjmmWeM/ppE1oLdc0REREQG4Og5IiIiIgMwaCIiIiIyAIMmIiIiIgMwaCIiIiIyAIMmIiIiIgMwaCIiIiIyAIMmIiIiIgMwaCIiIiIyAIMmIiIiIgMwaCIiIiIyAIMmIiIiItze/wNQOcS0TMIjcAAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import scipy.stats\n", "try:\n", " import matplotlib.pyplot as plt\n", "except ImportError:\n", " plt = None\n", " print('matplotlib is not installed: the plot is skipped.')\n", "\n", "if plt is not None:\n", " k = validation.decisions.index('lambda=0.9')\n", " x = validation.impacts[k]\n", " grid = np.linspace(x.min(), np.quantile(x, 0.9995), 400)\n", " fig, ax = plt.subplots(figsize=(6, 3.5))\n", " ax.hist(x, bins=200, range=(grid[0], grid[-1]), density=True, alpha=0.6, label='draws')\n", " ax.plot(grid, scipy.stats.norm.pdf(grid, point.mean, point.sigma), label='normal, closed-form moments')\n", " ax.axvline(point.adjusted, color='k', ls='--', label=r'$z^*$ at $\\lambda = 0.9$')\n", " ax.set_xlabel('impact X (kg CO$_2$-eq)'); ax.legend(); fig.tight_layout()" ] }, { "cell_type": "markdown", "id": "471cbb40", "metadata": {}, "source": [ "**When $X$ is exactly normal, the coverage is exact.** Replace every declared $B$\n", "entry by a normal with the same mean and variance and make the CFs exact: then $X$\n", "is normal, and the empirical coverage of the impact row must equal $\\lambda_z$ up to\n", "Monte Carlo error. This is the check that separates the method from its one\n", "approximation." ] }, { "cell_type": "code", "execution_count": 16, "id": "fa5f3e56", "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", " \n", " \n", " \n", " \n", " \n", "
lambda_impactcoveragesecoverage_lowcoverage_highskew
lambda
0.500.7500.75110.00100.74920.7530-0.0030
0.900.9500.95060.00050.94960.9515-0.0032
0.990.9950.99490.00020.99460.9952-0.0034
\n", "
" ], "text/plain": [ " lambda_impact coverage se coverage_low coverage_high skew\n", "lambda \n", "0.50 0.750 0.7511 0.0010 0.7492 0.7530 -0.0030\n", "0.90 0.950 0.9506 0.0005 0.9496 0.9515 -0.0032\n", "0.99 0.995 0.9949 0.0002 0.9946 0.9952 -0.0034" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import copy\n", "normal = copy.deepcopy(data)\n", "for db, block in normal['If'].items():\n", " specs = {}\n", " for index, spec in block['declared'].items():\n", " mean, var = unc.spec_moments(spec)\n", " specs[index] = {'uncertainty_type': 3, 'loc': mean, 'scale': var ** 0.5}\n", " unc.override(normal, 'If', db, specs)\n", "unc.override(normal, 'Cf', method, {e: {'uncertainty_type': 1} for e in normal['Cf'][method]['declared']})\n", "\n", "normal_problem = unc.ChanceConstrained(worker, unc.compute_moments(normal, worker),\n", " upper_bounds={electrolysis: capacity})\n", "normal_front = normal_problem.solve([0.5, 0.9, 0.99])\n", "check = unc.validate(normal_front, normal_problem, normal, n=200_000, seed=7).table.loc['front']\n", "check.assign(se=np.sqrt(check['lambda_impact'] * (1 - check['lambda_impact']) / 200_000))[\n", " ['lambda_impact', 'coverage', 'se', 'coverage_low', 'coverage_high', 'skew']]" ] }, { "cell_type": "markdown", "id": "2b6c39fd", "metadata": {}, "source": [ "## 7. Building blocks\n", "\n", "The reduced system and the sampler are public, so a study can formulate its own\n", "problem on top of them, e.g. a scenario-based CVaR front. `problem.projections()`\n", "returns $s_0$, $\\mathrm{E}[X]=m_0+m^\\top v$, the rows of $S$ on the processes that\n", "carry variance and $\\mathrm{E}[B_u]S$; `draw_parameters` returns draws of the declared\n", "parameters as arrays. Together they give the sampled impact of the base and of\n", "one unit of each alternative without further solves." ] }, { "cell_type": "code", "execution_count": 17, "id": "68f779bd", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "m0 = 1.8659, m = [10.634 8.8386], S_J (4, 2), B_u S (1, 2), draws: B (7, 1000), CFs (2, 1000)\n" ] } ], "source": [ "p = problem.projections()\n", "draws = unc.draw_parameters(data, 1000, seed=1, processes=p.J)\n", "print(f\"m0 = {p.m0:.4f}, m = {np.round(p.m, 4)}, S_J {p.S_J.shape}, B_u S {p.B_unc_S.shape}, \"\n", " f\"draws: B {draws.B.shape}, CFs {draws.Q.shape}\")" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "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.12.13" } }, "nbformat": 4, "nbformat_minor": 5 }