{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "84809625",
   "metadata": {},
   "source": [
    "# Ex 1: Analytical Reference Model (2D / 4-DOF)\n",
    "\n",
    "This notebook establishes the baseline statistical tolerance analysis for the minimal 1.5D/2D assembly: a rigid square part that must fit within a square pocket. \n",
    "\n",
    "**Hypotheses & Setup (from Thesis Section 4.2.1):**\n",
    "- **Dimensions:** Male part $L_1 = 99.8$ mm, Female pocket $L_2 = 100$ mm, Feature width $h = 10$ mm.\n",
    "- **Tolerance:** $t = 0.31$ mm applied to the two interacting features.\n",
    "- **Capability:** Target process capability $C_p = 1$.\n",
    "- **Degrees of Freedom:** 4 variables in total. \n",
    "  - $u_{d4}, \\gamma_{d4}$ for the Planar Feature on Part 2 (Pocket).\n",
    "  - $u_{d5}, \\gamma_{d5}$ for the Planar Feature on Part 1 (Square Part).\n",
    "- **Limit State:** The standard physical clearance $j$ must remain positive. The continuous slack variable transformation ($s$) is not required for this purely analytical Monte Carlo evaluation.\n",
    "\n",
    "\n",
    "\n",
    "<div>\n",
    "<img src=\"../Pictures/Ex1_Rectangular_Pocket.png\" width=\"750\"/>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "922dfb66-e8d7-4fa2-8b96-4508afd9dd0a",
   "metadata": {},
   "source": [
    "## Defining the \"Max Sigma\" Baseline\n",
    "\n",
    "Before introducing epistemic uncertainty and exploring the credal set, we must establish a precise probabilistic baseline. We evaluate the assembly at the boundary of the capability constraint. This assumes the manufacturing process uses the maximum allowable variance dictated by $t$ and $Cq_p$.\n",
    "\n",
    "For a planar feature, the maximal standard deviations for translation and rotation, assuming independent centered normal distributions, are extracted from the local defect expression:\n",
    "$$ \\sigma_{u, \\max} = \\frac{t}{6 C_p} $$\n",
    "$$ \\sigma_{\\gamma, \\max} = \\frac{t}{6 C_p \\cdot y_{\\max}} $$\n",
    "where $y_{\\max} = h/2$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "bc12c066",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Max Translational Std Dev (sigma_u_max): 0.05167\n",
      "Max Rotational Std Dev (sigma_gamma_max): 0.01033\n"
     ]
    }
   ],
   "source": [
    "import openturns as ot\n",
    "import numpy as np\n",
    "import otaf\n",
    "\n",
    "# Geometric parameters\n",
    "L1 = 99.8  # Nominal length of the male piece (Part 1)\n",
    "L2 = 100.0 # Nominal length of the female pocket (Part 2)\n",
    "h = 10.0   # Nominal width of the features\n",
    "y_max = h / 2.0\n",
    "\n",
    "# Tolerance and Capability\n",
    "t = 0.31\n",
    "Cqp = 1.0\n",
    "\n",
    "# Maximal allowable standard deviations based on the capability constraint\n",
    "sigma_target = t / (6 * Cqp)\n",
    "\n",
    "sigma_u_max = sigma_target\n",
    "sigma_gamma_max = sigma_target / y_max\n",
    "\n",
    "print(f\"Max Translational Std Dev (sigma_u_max): {sigma_u_max:.5f}\")\n",
    "print(f\"Max Rotational Std Dev (sigma_gamma_max): {sigma_gamma_max:.5f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "99532646-3b1d-430d-aa71-2d9148f2ad73",
   "metadata": {},
   "source": [
    "## Analytical Assembly Function\n",
    "\n",
    "The assembly constraint is the clearance $j = X_2 - X_1$, which must remain strictly positive across the entire width of the feature. Because the parts can rotate, the minimum clearance will always occur at one of the extremities ($+y_{\\max}$ or $-y_{\\max}$). \n",
    "\n",
    "The function calculates the physical position of the top and bottom corners of both features and extracts the minimum resulting gap."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "9ddf68f7-10b5-4884-a973-fd44b629b5aa",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def analytical_assembly_model_4_DOF(sample_of_defects, round_val=12):\n",
    "    \"\"\"\n",
    "    Calculate the minimum spatial clearance between the two parts.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    sample_of_defects : ot.Sample\n",
    "        OpenTURNS sample containing [u_d4, gamma_d4, u_d5, gamma_d5]\n",
    "    round_val : int\n",
    "        Decimal rounding to handle floating point noise.\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    ot.Sample\n",
    "        The minimum gap j for each realization.\n",
    "    \"\"\"\n",
    "    size = sample_of_defects.getSize()\n",
    "    desc = sample_of_defects.getDescription()\n",
    "    \n",
    "    # Safely extract variables based on thesis nomenclature\n",
    "    zero_array = np.zeros((size, 1))\n",
    "    \n",
    "    smp_u4     = np.array(sample_of_defects.getMarginal([\"u_d4\"]))     if \"u_d4\" in desc else zero_array\n",
    "    smp_gamma4 = np.array(sample_of_defects.getMarginal([\"gamma_d4\"])) if \"gamma_d4\" in desc else zero_array\n",
    "    \n",
    "    smp_u5     = np.array(sample_of_defects.getMarginal([\"u_d5\"]))     if \"u_d5\" in desc else zero_array\n",
    "    smp_gamma5 = np.array(sample_of_defects.getMarginal([\"gamma_d5\"])) if \"gamma_d5\" in desc else zero_array\n",
    "\n",
    "    # Part 1 (Male) displaced geometry at extremities\n",
    "    # Displacement = nominal + translation +/- lever arm effect of rotation\n",
    "    X1_corners = np.asarray([\n",
    "        L1 + smp_u5 - y_max * smp_gamma5, \n",
    "        L1 + smp_u5 + y_max * smp_gamma5\n",
    "    ]) \n",
    "\n",
    "    # Part 2 (Female) displaced geometry at extremities\n",
    "    X2_corners = np.asarray([\n",
    "        L2 - smp_u4 - y_max * smp_gamma4, \n",
    "        L2 - smp_u4 + y_max * smp_gamma4\n",
    "    ])\n",
    "    \n",
    "    # The physical gap is evaluated at both corners, the minimum dictates collision\n",
    "    gaps = X2_corners - X1_corners\n",
    "    min_gap = np.expand_dims(np.squeeze(gaps.min(axis=0)), axis=1).round(round_val).tolist()\n",
    "    \n",
    "    return ot.Sample(min_gap)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "173208cb-79ca-47fd-b83c-6a2b21048c1b",
   "metadata": {},
   "source": [
    "## Monte Carlo Propagation\n",
    "\n",
    "We estimate the reference probability of failure. According to the equal allocation hypothesis, this maximum standard deviation configuration should yield a $P_f \\approx 5\\%$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "c17de2b0",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING:root:No mu_dict passed, initializing all unspecified means to 0.0\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Reference Probability of Failure at Max Sigma: 5.1710%\n"
     ]
    }
   ],
   "source": [
    "# Fixed seed for reproducible reference\n",
    "ot.RandomGenerator.SetSeed(888)\n",
    "size_MC = int(1e5)\n",
    "\n",
    "# Constructing the distribution using the otaf library and thesis nomenclature\n",
    "deviation_vector = otaf.distribution.get_composed_normal_defect_distribution(\n",
    "    defect_names=[\"u_d4\", \"gamma_d4\", \"u_d5\", \"gamma_d5\"],\n",
    "    sigma_dict={\"u_\": sigma_u_max, \"gamma_\": sigma_gamma_max}\n",
    ")\n",
    "\n",
    "# Sample and evaluate\n",
    "max_defect_sample = deviation_vector.getSample(size_MC)\n",
    "clearance_sample = analytical_assembly_model_4_DOF(max_defect_sample)\n",
    "clearance_distribution = ot.UserDefined(clearance_sample)\n",
    "\n",
    "# Compute Probability of Failure (Gap <= 0)\n",
    "pf_baseline = clearance_distribution.computeCDF(0.0)\n",
    "print(f\"Reference Probability of Failure at Max Sigma: {pf_baseline:.4%}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "788cbc89-ca88-46a4-996c-ac8ca48f56fa",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAoAAAAHgCAYAAAA10dzkAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAATRZJREFUeJzt3X98zfX///H72W8z86OxGfOjhPSDEK13indDJeXdD8Ibb5VKKayEkpHepvKrHyTK2/v9jpBKfb4kS6ZEFFH5VYq3wsbyY2Ozne28vn+82mZtZ2y289o5r9v1cjmXPV+v83qd8ziP5nTf66fDMAxDAAAAsA0/qwsAAACAZxEAAQAAbIYACAAAYDMEQAAAAJshAAIAANgMARAAAMBmCIAAAAA2QwAEAACwGQIgAACAzRAAAQAAbIYACMDrnDhxQg8++KDq1q2r6tWrq0uXLtq6dWuZXuODDz7QLbfcooiICAUFBSk6Olq9e/fWZ599VrBMcnKyHA5HwSM4OFiRkZHq3LmzJk+erKNHjxZ73QULFhRZ5+zHmDFjLvizA0BFCLC6AAAoC5fLpR49emj79u0aNWqUIiIiNHv2bHXu3FlbtmzRpZdeWur6hmHovvvu04IFC3T11VcrPj5eUVFROnz4sD744APddNNN+vLLL3XdddcVrPP444/rmmuuUV5eno4ePaoNGzYoISFB06dP19KlS/XXv/612Ps899xzatq0aZF5V1xxRcU0AQAuEAEQgFdZtmyZNmzYoHfffVd33323JKl3795q3ry5EhIStGjRolLXnzZtmhYsWKARI0Zo+vTpcjgcBc8988wz+u9//6uAgKJfjZ06dSp4r3zbt29Xt27ddNddd2nnzp2qX79+kedvueUWtW/f/kI+KgBUGnYBA7DcsmXL5HA4tG7dumLPvfHGG3I4HPrhhx8Klo2MjNSdd95ZsEzdunXVu3dvffjhh8rOznb7PllZWUpMTFTLli01derUIuEv34ABA9ShQ4dz1ty6dWvNnDlTJ06c0GuvvXY+HxMAqgwCIADL9ejRQ2FhYVq6dGmx55YsWaLLL7+8YPfpt99+q7Zt28rPr+jXV4cOHZSZmakff/zR7fusX79ex44dU79+/eTv73/Bdd99992qVq2aVq9eXey5kydPKi0trcgDAKoKAiAAy1WrVk09e/bUsmXLlJeXVzA/JSVF69atU58+fQrmHT58uNjuVkkF8w4dOuT2fXbt2iVJuvLKKyuk7sDAQDVv3lw///xzsefi4uJUt27dIg8AqCo4BhBAldCnTx+98847Sk5O1k033STJ3N3rcrmKBMCsrCwFBwcXWz8kJKTgeXfS09MlSTVq1KiwusPCwpSRkVFs/qxZs9S8efMKex8AqEgEQABVws0336yaNWtqyZIlBQFwyZIlatOmTZEgVa1atRKP8ztz5kzB8+6Eh4dLUomBrbxOnTpVYqDs0KEDJ4EAqLLYBQygSggODlavXr30wQcfKDc3VwcPHtSXX35ZZOufZO7qPXz4cLH18+dFR0e7fY+WLVtKkr7//vsKqdnpdOrHH39Us2bNKuT1AMBTCIAAqow+ffooLS1Na9as0bvvvivDMIoFwDZt2mjr1q1yuVxF5m/atEmhoaGl7na9/vrrVbt2bb3zzjtFjjUsr2XLlikrK0vdu3e/4NcCAE8iAAKoMuLi4lSnTh0tWbJES5YsUYcOHYpdTPnuu+9Wamqq3n///YJ5aWlpevfdd9WzZ88Sjw/MFxoaqtGjR2vXrl0aPXq0DMMotszbb7+tzZs3n7PW7du3a8SIEapdu7YeffTRMnxKALAexwACqDICAwN15513avHixTp9+rSmTp1abJm7775b1157rQYPHqydO3cW3AkkLy9PEydOPOd7jBo1Sjt27NC0adO0du1a3X333YqKilJKSoqWL1+uzZs3a8OGDUXW+eKLL3TmzBnl5eXp999/15dffqmPPvpINWvW1AcffKCoqKgK6wEAeILDKOlPYACwyKeffqquXbvK4XDowIEDatiwYbFljh8/rlGjRmn58uXKysrSNddco6lTp5bppIv33ntPc+fO1TfffKP09HTVrVtXN9xwg4YOHaobb7xRknkv4C5duhSsExgYqFq1aumyyy5T9+7dNWTIkGKXd1mwYIEGDx6sr7/+mpNAAFRZBEAAAACb4RhAAAAAmyEAAgAA2AwBEAAAwGYIgAAAADZDAAQAALAZAiAAAIDNEAABAABshgAIAABgMwRAAAAAmyEAAgAA2AwBEAAAwGYIgAAAADZDAAQAALAZAiAAAIDNEAABAABshgAIAABgMwRAAAAAmyEAAgAA2AwBEAAAwGYIgAAAADZDAAQAALAZAiAAAIDNEAABAABshgAIAABgMwRAAAAAmyEAAgAA2AwBEAAAwGYIgAAAADZDAAQAALAZAiAAAIDNEAABAABshgAIAABgMwFWF+DNXC6XDh06pBo1asjhcFhdDgAAOA+GYSgjI0PR0dHy87PntjAC4AU4dOiQYmJirC4DAACUw6+//qqGDRtaXYYlCIAXoEaNGpLMX6Dw8HBJktPp1OrVq9WtWzcFBgZaWV6VQ29KR3/cozfu0Rv36I17du9Nenq6YmJiCv4/bkcEwAuQv9s3PDy8SAAMDQ1VeHi4Lf9RlYbelI7+uEdv3KM37tEb9+iNyc6Hb9lzxzcAAICNEQABAABshgAIAABgMxwDWMkMw5DT6VRubq7VpVjO6XQqMDBQmZmZHj3mJCAgQIGBgbY+1gMAgLMRACtRdna29u/fr1OnTlldSpURGRmpvXv3evx9w8LC1KRJEwUHB3v8vQEAqGoIgJXE5XJp586dCggIUNOmTRUcHMwWKAsYhqHs7Gz99ttv2rFjh1q2bKnQ0FCrywIAwFIEwEqSnZ0tl8ulpk2bKiwszOpybK169eoKCgrSnj17tGrVKnXt2tXW134CAICTQCqJYRiSZNtbzFQ1+f8dDh48qJUrV8rpdFpcEQAA1iGdwFbq1KmjQ4cO6eTJk1aXAgCAZXwmAH7++efq2bOnoqOj5XA4tHz58nOuk5ycrLZt2yo4OFjNmjXTggULKr1OWCsgIEB5eXnKzs62uhQAACzjMwHw9OnTat26tWbNmnVey+/bt089evRQly5dtG3bNo0YMUIPPPCAPvnkk0qu1PuVJzgbhqG5c+eqY8eOCgsLU61atdS+fXvNnDlTmZmZkqQJEybI4XDI4XAoICBAERERuuGGGzRz5sxiga1z584Fy579OJ/L7eTvngcAwK585iSQW265Rbfccst5Lz9nzhw1bdpU06ZNkyRddtllWr9+vWbMmKHu3btXVpleLz84P/zww1q4cKHWrFmjBx54QPXr1y+1bwMGDND777+vcePG6bXXXlPdunW1fft2zZw5U02aNFGvXr0kSZdffrk+/fRTuVwu/f7770pOTtbzzz+v//73v0pOTi5y8saQIUP03HPPFXmfgACf+ZUGAKDS2Pb/lhs3blRcXFyRed27d9eIESPcrpOdnV1kS1R6erok8wLH+ScV5P/0xgs/z507VxMmTNBvv/1W5OSVO+64QxdddJHmz59fruC8dOlSLVy4UMuXL9cdd9xRML9Jkya6/fbbC/oomQEuKipKkhQdHa0rr7xSXbt2VevWrfXCCy/o+eefL1g2NDS0YNnzZRiGXC6XcnNzq9yJIH/+HUIheuMevXHPa3vjckmnT0snT0pZWVJmphyZmea848el9HTJz0+OP55TVpaUkyPHb79JGRky6tWTw+mU/ng4vv9eRni4dNFFUl6elJcnP6dTf0lLk9+UKXK5XPLbvFlGdLRUp07BMnK5Ch9nzXOkpJiv5+9v1msY5uNc47Ise/a4RQvlfvddhbbY634nKoFtA2BKSooiIyOLzIuMjFR6erqysrJUrVq1YuskJiZq4sSJxeavXr262LXlNm3aVOz1q7p77rlHjz32mNauXaubbrpJknTs2DGtWrVKK1eulFS+4Lxw4UK1aNGiSPjL53A4VLNmzVLratmypW655Ra9//77RQJgefz66686fPiw1q1bd873tUpSUpLVJVRZ9MY9euOelb1xOJ0KTk9X8PHjCvnjEXz8uIIyMhR06pQCT51SUEaG6uzZo5ywMAVV0o0DSroKbcSflzl0SDp06Pxe76w/3CtbRkaGPvvj/0EVJf/QIzuzbQAsj7Fjxyo+Pr5gOj09XTExMerWrZvCw8MlmX9VJCUlqWPHjtq/f3/xF2nfXkpJ8VDFf4iKkr755pyL1a5dW7fccosWLVpUEACXLVumiIgIdenSRVL5gvNPP/2kFi1aXNBHaNmypVavXl1k3uzZs/Xmm28WTD/00EMFWybdiYmJkcvl0o033qjo6OgLqqmi5f/udO3a1aO3yvMG9MY9euOeR3qTni79/LMcv/wix9690uHDchw4YG6NO3xYjtTU836pygp/58Pw95cjL88ch4WZW/f8/MyfJYwd+/aZy156qfkCDof5ONf4HM8bJcyv3rixbr311gr9vOkeDLBVlW0DYFRUlFL/9A8zNTVV4eHhJYYYSQoODi7xVmKBgYHFvlzcHouWkiIdPFi+oj2gf//+GjJkiGbPnq3g4GAtXLhQ99577wVdz7AiTrowDKPYnVT69++vZ555pmC6Vq1a53wdh8MhPz+/gvsDV0Ul/T7BRG/cozfuXXBvnE7p+++lrVul3bulffvMx6+/SmlpFVNk9epSbq6UnS1deql09Kh07bXS//4n3XCDFBpqLlOtmnTmjBQeLkVGmvOrVZOCg82HJIWFSUFBUmBg4cPfv/BnQICcLpc+Xr1at9x2mwL/dKeqstyzqqLvb+Xu9Sr6jFX+rdg4AMbGxhbs1syXlJSk2NjYyn3jMh6z5un37NmzpwzD0IoVK3TNNdfoiy++0IwZM856qbIH5+bNm2v37t3lq/0Pu3btUtOmTYvMq1mzppo1a3ZBrwsARRiGGfK+/NLcc/L119LOnWboKgt/f/O7t2FD82f+o35981G3rlS7tnnMXa1aZmDzJKdTRn4g5DaltuQzAfDUqVPau3dvwfS+ffu0bds21alTR40aNdLYsWN18OBB/ec//5EkPfzww3rttdf01FNP6b777tNnn32mpUuXasWKFZVb6HnsirVSSEiI7rzzTi1cuFB79+5VixYt1LZt24LnyxOc+/Xrp3vvvVcffvhhseMADcNQenp6qcfj7d69W6tWrdLYsWPL+akAwI2sLGnTJmn9emndOvM7+sSJc6/n7y81aCA1amRusbvkEvNn48ZSdLT5yD9JAqiCfCYAfvPNNwXHqUkqOFZv0KBBWrBggQ4fPqwDBw4UPN+0aVOtWLFCI0eO1Msvv6yGDRvqzTff5BIwMnet3nbbbdqxY4f+/ve/F3muPMG5d+/e+uCDD9S3b1+NGzdO3bp1U926dfX9999rxowZeuyxxwouA5Obm6uUlJRil4Fp06aNRo0aVZkfG4BdHD0qrVghvfee9OmnpW/d8/OTmjWT2rUzH23aSBdfbG7ZYzcivJjPBMDOnTuXeqxZSRcr7ty5s7799ttKrMo7/fWvf1WdOnW0Z88e9evXr8hz5QnODodDixYt0ty5czV//nz985//VEBAgC699FINHDiwyLo7duxQ/fr15e/vr5o1a6pVq1YaO3ashg4dWuLxlwBwXo4dkxYtkpYvl5KTzUualCQqygx6N94odehgjsPCPFkp4BE+EwBRcfz8/HSolEsBlCc4+/n56eGHH9bDDz/sdpkJEyZowoQJ5/V6ycnJZXp/ADZ0/LiarFol/zlzzC19JV2fNTJS6tFD6tTJDH1NmnBMHGyBAAgA8C0//iglJirgnXfUuqT7fjdtKt19t9SrlxQbS+CDLREAAQDezzCkDRukadPM3byGUfSSItHRUu/e0uDB0pVXEvpgewRAAID3cjqlZcvM4LdlS5GnjJo19UunTmr07LMKvOYaQh9wloq+tiIAAJUvJ0dasEC67DKpX7+i4a9+fSkxUbm//KIfHnhAuvpqwh/wJ2wBBAB4j8xMae5caerU4ndVattWevxx6d57zbtiOJ3W1Ah4AQJgJcm/rY7L5bK4Ekj8dwC8Xna29MYb0pQp0uHDRZ/r0kV65hnpr39lSx9wngiAlST/PoOnTp1SGNeQstypP26ynpOTY3ElAMrsgw+kYcOkP1+e6vbbpdGjpeuus6YuwIsRACtJQECAIiIidPCPXRRhYWHy8+OQS09zuVw6deqUDh48qBMnTsjlcsnhcBS58TmAKurgQWnoUOn//q/o/DvvlMaPl1q3tqYuwAcQACtRo0aNJKkgBMI6J06cUGpqqjIzMxUYGKhq1apZXRIAd3Jzpddfl559Vjp5snB+165SYqJ5dw4AF4QAWIkcDocaN26szMxMffnll8rNzVX16tWtLssyhmHo0KFDio6O9ugWuJycHOXm5srpdCojI0Pt27dXeHi4x94fQBls3So98IB09t2GoqKkGTOkPn04xg+oIARAD7jssstkGIY2bNigkydP2vaEBMMwlJaWpmrVqnl8F6zD4VBAQICuueYade7cWf7+/h59fwDnkJcnTZpkPs7+jhw0SJo+XapTx7raAB9EAPSQVq1aqWXLlsrKypJhGFaXYwmn06mkpCR17dq14CQZTwoODrbkfQGcwy+/mEFv/frCeVdeKc2aZd6jF0CFIwB6kJ+fn613ATudTgUHByssLIwgBsC0YoX0979LJ06Y035+5gkeY8dKQUGWlgb4MgIgAMDzsrOlESOkOXMK5zVqJP3731LnzlZVBdgGARAA4FlpadI990jJyYXzevUywx8naAEewYXpAACe89VX0rXXFoa/kBDzWL/33yf8AR7EFkAAgGcsWWKe7JGdbU5HRZnBLzbW2roAG2ILIACg8r31ltS3b2H469hR2riR8AdYhAAIAKg8hiE984x5cef8S2Ddd5/0+edSkyaWlgbYGbuAAQCVIzdXGjZMeuONwnkjRpgXduaOHoClCIAAgIqXk2Peum358sJ5L78sPfYY4Q+oAgiAAICKlZUl3XmntGqVOR0YKM2fb17wGUCVQAAEAFScEyekHj2kDRvM6WrVzK2A3bpZWRWAPyEAAgAqxrFjUvfu0jffmNNhYeat3m64wdq6ABRDAAQAXLgjR6S4OOn7783piy6SVq+W2ra1ti4AJSIAAgAuTHq61LNnYfiLjJSSkqQrr7S2LgBuEQABAOWXlWXex3fzZnM6Olpau1Zq3tzSsgCUjgtBAwDKJyfHPNt37VpzunZt6eOPCX+AFyAAAgDKZ9iwwku9hIVJ//d/0lVXWVsTgPNCAAQAlN3kydK8eeY4JET66CPpL3+xtiYA540ACAAom7ffNu/vm++NN6QuXayrB0CZEQABAOdv2zbpvvsKpxMTpYEDLSsHQPkQAAEA5ycjQ+rfX3I6zemHHpJGj7a2JgDlQgAEAJxbXp7Ut6+0c6c53aaN9MorksNhaVkAyocACAAonWFIjz9u3tZNkmrVkt55RwoKsrQsAOVHAAQAlG7uXGn2bHPs7y+9+67UsqW1NQG4IARAAIB7n35qXu8v31tvmff8BeDVCIAAgJL9/LPUp4+Um2tOjxghDRpkaUkAKgYBEABQ3JkzUu/e0rFj5nSPHtLUqdbWBKDCEAABAMUNHy5t3WqOL73UvPizv7+1NQGoMARAAEBRy5aZJ35IUrVq5kkftWpZWhKAikUABAAU2rVLGjy4cPrll6XWra2rB0ClIAACAEwZGdLdd0unTpnT994rPfCAtTUBqBQEQACA6fHHC+/0cfnl5iVfuNMH4JMIgAAAaelSacECcxwWJr33nhQaamlJACoPARAA7O6336SHHy6cfvVVqUUL6+oBUOkIgABgZ4YhPfSQdPy4Od27Nxd7BmyAAAgAdjZ/vrRypTmuX196/XWO+wNsgAAIAHb1v/9JTzxROD17tlSnjnX1APAYAiAA2JHLJT34oHTypDk9YIDUq5elJQHwHAIgANjR3LnS6tXmODraPPEDgG0QAAHAbg4flsaOLZx+6y2pZk3r6gHgcQRAALCb8eOlEyfM8YAB0s03W1oOAM8jAAKAnWzZYp75K0nh4dK0adbWA8ASBEAAsIvcXPOafy6XOf3001LdutbWBMASPhUAZ82apSZNmigkJEQdO3bU5s2bS11+5syZatGihapVq6aYmBiNHDlSZ86c8VC1AOBhL79sbgGUzHv9xsdbWw8Ay/hMAFyyZIni4+OVkJCgrVu3qnXr1urevbuOHDlS4vKLFi3SmDFjlJCQoF27dumtt97SkiVL9PTTT3u4cgDwgP37pWefLZyeN08KDLSsHADW8pkAOH36dA0ZMkSDBw9Wq1atNGfOHIWGhmp+/rEuf7Jhwwb95S9/Ub9+/dSkSRN169ZNffv2PedWQwDwSg8+KGVlmeNHHpFiY62tB4ClfCIA5uTkaMuWLYqLiyuY5+fnp7i4OG3cuLHEda677jpt2bKlIPD98ssvWrlypW699VaP1AwAHvPhh1JSkjmOiZH++U9r6wFguQCrC6gIaWlpysvLU2RkZJH5kZGR2r17d4nr9OvXT2lpabr++utlGIZyc3P18MMPl7oLODs7W9nZ2QXT6enpkiSn0ymn01kwPvsnCtGb0tEf9+iNe+fszenTCnjsMeXf3Tc3MVFG9eqSDXrJ7417du+NXT/32XwiAJZHcnKyJk+erNmzZ6tjx47au3evhg8frkmTJunZs4+TOUtiYqImTpxYbP7q1asVGhpaZF5S/l/bKIbelI7+uEdv3HPXmxbvvKOWv/4qSTrSurU2Vq8urVzpydIsx++Ne3btTWZmptUlWM5hGIZhdREXKicnR6GhoVq2bJl6nXUvy0GDBunEiRP68MMPi63TqVMnXXvttXrppZcK5r399tt68MEHderUKfn5Fd87XtIWwJiYGKWlpSk8PFyS+VdFUlKSunbtqkAOsC6C3pSO/rhHb9wrtTd79yqgTRs5cnJkBAQod9s2qXlzS+q0Ar837tm9N+np6YqIiNDJkycL/v9tNz6xBTAoKEjt2rXTmjVrCgKgy+XSmjVrNGzYsBLXyczMLBby/P39JUnuMnFwcLCCg4OLzQ8MDCz2D6ikeTDRm9LRH/fojXsl9mb8eCknR5LkGDFCgZdfbkFl1uP3xj279saOn/nPfCIASlJ8fLwGDRqk9u3bq0OHDpo5c6ZOnz6twYMHS5IGDhyoBg0aKDExUZLUs2dPTZ8+XVdffXXBLuBnn31WPXv2LAiCAOC1Pv1UWrbMHEdEmGEQAP7gMwGwT58+Onr0qMaPH6+UlBS1adNGq1atKjgx5MCBA0W2+I0bN04Oh0Pjxo3TwYMHVbduXfXs2VP/5Ow4AN7O6ZSeeKJw+qWXpBo1rKsHQJXjMwFQkoYNG+Z2l29ycnKR6YCAACUkJCghIcEDlQGAB73xhvTdd+b46qulgQOtrQdAleMT1wEEAPzh5Mmiu3tnzZJKOKkNgL3xrQAAviQxUTp+3Bz3788dPwCUiAAIAL7i11+ll182x0FB0vPPW1sPgCqLAAgAvmLiROnMGXP8+ONSkyaWlgOg6iIAAoAv2L5dmj/fHIeHS089ZW09AKo0AiAAeDvDkP/IkVL+RezHjZPq1rW2JgBVGgEQALxcw3Xr5Ld+vTlxySXSY49ZWxCAKo8ACADe7MQJXf6vfxVOv/KKFBJiXT0AvIJPXQgaAOzGb/JkBZ48aU706iXdequl9QDwDmwBBABvtWuX/F57TZJkBAcXXgIGAM6BAAgA3io+Xo7cXEmS68knpUaNLC4IgLcgAAKAN1q9Wlq1SpKUWbeuXFz2BUAZEAABwNvk5UkjRxZM7u7XT6pWzcKCAHgbAiAAeJsFC6SdOyVJrg4d9OuNN1pbDwCvQwAEAG+SmSmNH18w6ZoyRfLjqxxA2fCtAQDe5NVXpUOHzPFtt8m4/npr6wHglQiAAOAt0tOladPMsZ+f9MIL1tYDwGsRAAHAW8ycKR09ao5795ZatbK0HADeiwAIAN4gNdUMgJK59e/55y0tB4B3IwACgDeYPl06ftwc//3v0iWXWFsPAK9GAASAqu6336RXXjHHQUFSYqK19QDwegRAAKjqxo2Tzpwxx8OGSdHR1tYDwOsRAAGgKvvpJ+m//zXHtWpJTz9taTkAfAMBEACqsvHjJZfLHD/5pHTRRdbWA8AnEAABoKr67jtpyRJzXLeu9Pjj1tYDwGcQAAGgqvrnPyXDMMdjxkg1alhbDwCfQQAEgKrou++kd981x/XqSUOHWlsPAJ9CAASAqmjMmMKtf6NGSdWqWVsPAJ9CAASAqmbzZunjj81xo0bmpV8AoAIRAAGgqhk/vnD89NNSSIh1tQDwSQRAAKhKkpOlTz4xx40bS4MHW1oOAN9EAASAquS554qOg4KsqwWAzyIAAkBVsXat+ZCkZs2kfv2srQeAzyIAAkBVcfaxfwkJUkCAdbUA8GkEQACoCjZvltavN8ctW0p9+1pbDwCfRgAEgKogMbFw/MQTkr+/dbUA8HkEQACw2nffScuXm+P69aUBAywtB4DvIwACgNXO3vo3dqwUHGxdLQBsgQAIAFbav7/wnr9160r3329pOQDsgQAIAFaaMUPKyzPHjz4qhYZaWw8AWyAAAoBVjh6V5s41x9WqmQEQADyAAAgAVnnlFenMGXP88MNSRIS19QCwDQIgAFghM1N6/XVzHBgojRhhaTkA7IUACABWmD9f+v13c3zPPVKjRtbWA8BWCIAA4Gk5OdJLLxVOP/WUdbUAsCUCIAB42n/+Ix04YI5vvVVq3draegDYDgEQADzJ5ZKmTSucfvZZ62oBYFsEQADwpKVLpd27zXGnTtK111pbDwBbIgACgKcYhjR5cuH0uHHW1QLA1giAAOApa9ZI339vjjt2lLp2tbYeALZFAAQAT5kypXAcHy85HNbVAsDWCIAA4AmbNplbACXpkkuku+6yth4AtkYABABPmDixcDxmjOTvb10tAGyPAAgAle3HH6WPPzbHjRtLAwdaWw8A2yMAAkBlmzq1cDxsmBQUZF0tACACIABUrv/9T1qwwByHh0v3329pOQAg+VgAnDVrlpo0aaKQkBB17NhRmzdvLnX5EydO6NFHH1X9+vUVHBys5s2ba+XKlR6qFoAtvPSS5HSa48cek2rXtrYeAJAUYHUBFWXJkiWKj4/XnDlz1LFjR82cOVPdu3fXnj17VK9evWLL5+TkqGvXrqpXr56WLVumBg0a6H//+59q1arl+eIB+Ka0NOlf/zLHoaHmpV8AoArwmQA4ffp0DRkyRIMHD5YkzZkzRytWrND8+fM1ZsyYYsvPnz9fx44d04YNGxQYGChJatKkiSdLBuDrZs+WMjPN8X33SXXqWFsPAPzBJ3YB5+TkaMuWLYqLiyuY5+fnp7i4OG3cuLHEdT766CPFxsbq0UcfVWRkpK644gpNnjxZeXl5niobgC/LypJmzTLHfn7Sk09aWw8AnMUntgCmpaUpLy9PkZGRReZHRkZqd/5N1//kl19+0Weffab+/ftr5cqV2rt3rx555BE5nU4lJCSUuE52drays7MLptPT0yVJTqdTzj+O8fnzTxSiN6WjP+55Y2/83nxT/keOSJJcvXopLzq68FjACuSNvfEUeuOe3Xtj1899NodhGIbVRVyoQ4cOqUGDBtqwYYNiY2ML5j/11FNat26dNm3aVGyd5s2b68yZM9q3b5/8/7gg6/Tp0/XSSy/p8OHDJb7PhAkTNPHsi7n+YdGiRQoNDa2gTwPA2zny8nTT0KGq/kcATJ42TScvucTiqgDky8zMVL9+/XTy5EmFh4dbXY4lfGILYEREhPz9/ZWamlpkfmpqqqKiokpcp379+goMDCwIf5J02WWXKSUlRTk5OQoq4TpdY8eOVfxZB3Gnp6crJiZG3bp1K/gFcjqdSkpKUteuXQuOLYSJ3pSO/rjnbb1xLFumgPytf1276i+PPVZp7+VtvfEkeuOe3XuTvwfPznwiAAYFBaldu3Zas2aNevXqJUlyuVxas2aNhg0bVuI6f/nLX7Ro0SK5XC75+ZmHQv7444+qX79+ieFPkoKDgxUcHFxsfmBgYLF/QCXNg4nelI7+uOcVvTEMacaMgkm/UaPk54GavaI3FqE37tm1N3b8zH/mEyeBSFJ8fLzmzZunf//739q1a5eGDh2q06dPF5wVPHDgQI0dO7Zg+aFDh+rYsWMaPny4fvzxR61YsUKTJ0/Wo48+atVHAOAL1q2Tvv7aHLduLZ11choAVBU+sQVQkvr06aOjR49q/PjxSklJUZs2bbRq1aqCE0MOHDhQsKVPkmJiYvTJJ59o5MiRuuqqq9SgQQMNHz5co0ePtuojAPAFM2cWjp96SnI4LCsFANzxmQAoScOGDXO7yzc5ObnYvNjYWH311VeVXBUA2/jlF+n//s8cR0dL99xjbT0A4IbP7AIGAMvNmCG5XOb4kUckjjMCUEURAAGgIvz+u/Tmm+a4WjXpoYesrQcASkEABICKMG+edOaMOX7wQSkiwtp6AKAUBEAAuFB5edLrr5tjh0OqxOv+AUBFIAACwIV67z3pwAFzfPPNEnf9AFDFEQAB4EK4XNJzzxVODx9uXS0AcJ4IgABwIT77TNqxwxxfe63UrZu19QDAeSAAAsCFmDatcBwfz4WfAXgFAiAAlNf27dKqVea4cWPpb3+zth4AOE8EQAAorxdfLBw/+aQU4FM3VwLgwwiAAFAev/0mLVliji+6SLrvPmvrAYAyIAACQHm8/rp5/T9JGjpUCg21th4AKAMCIACUVXp64YWfAwLMAAgAXoQACABl9cor0vHj5rh/fyk62tp6AKCMCIAAUBanT0svv2yO/f2lZ56xth4AKAcCIACUxTvvSGlp5rhPH+nSS62tBwDKwZIAOHDgQGVkZBRMb9++XU6n04pSAOD8uVzS9OmF0yNGWFYKAFwISwLgwoULlZWVVTDdqVMn/frrr1aUAgDnLylJ2rXLHF9/vXTNNdbWAwDlZEkANAyj1GkAqJJmzCgcjxxpXR0AcIE4BhAAzseOHdInn5jjJk2k22+3tBwAuBCW3bdo586dSklJkWRuAdy9e7dOnTpVZJmrrrrKitIAoLipUwvHI0Zw2zcAXs2yb7CbbrqpyK7f2267TZLkcDhkGIYcDofy8q+yDwBW+u036e23zXGtWtz2DYDXsyQA7tu3z4q3BYDymTVLys01x48+KtWoYW09AHCBLAmAjRs3tuJtAaDsMjKkOXPMcWCgGQABwMtZehDLTz/9pA8//FD79++Xw+FQ06ZN1atXL1188cVWlgUAhRYskE6cMMf9+kn161tZDQBUCMsCYGJiosaPHy+Xy6V69erJMAwdPXpUY8aM0eTJk/Xkk09aVRoAmPLypJkzC6fj4y0rBQAqkiWXgVm7dq3GjRunZ555RmlpaTp8+LBSUlIKAuCYMWP0+eefW1EaABRasUL65RdzHBcncWUCAD7Cki2Ac+bM0QMPPKAJEyYUmV+nTh0999xzSklJ0euvv64bbrjBivIAwPTii4Xj4cOtqwMAKpglWwA3b96sAQMGuH1+wIAB+uqrrzxYEQD8yXffSV9+aY5btZJuvdXaegCgAlkSAFNTU9WkSRO3zzdt2rTgItEAYInXXiscP/KI5MeNkwD4Dku+0c6cOaOgoCC3zwcGBionJ8eDFQHAWdLSCi/8XKOGVMoeCwDwRpadBfzmm28qLCysxOcyMjI8XA0AnOW116SsLHM8eLAUHm5tPQBQwSwJgI0aNdK8efPOuQwAeFxGhvTKK+bY359LvwDwSZYEwP3791vxtgBwbm++KR0/bo779ZO4cxEAH2TJMYCfffaZWrVqpfT09GLPnTx5Updffrm++OILCyoDYGtnzkjTphVOjxljXS0AUIksCYAzZ87UkCFDFF7CcTU1a9bUQw89pOnTp1tQGQBbW7hQOnjQHN9+u3n5FwDwQZYEwO3bt+vmm292+3y3bt20ZcsWD1YEwPYMo+jWv7FjrasFACqZZdcBDAwMdPt8QECAjh496sGKANje2rXSrl3m+PrrpWuvtbYeAKhElgTABg0a6IcffnD7/Hfffaf69et7sCIAtvfCC4XjRx+1rg4A8ABLAuCtt96qZ599VmfOnCn2XFZWlhISEnTbbbdZUBkAW9qyRVq92hw3bSrddZe19QBAJbPkMjDjxo3T+++/r+bNm2vYsGFq0aKFJGn37t2aNWuW8vLy9Mwzz1hRGgA7eumlwvFTT0mlHKICAL7AkgAYGRmpDRs2aOjQoRo7dqwMw5AkORwOde/eXbNmzVJkZKQVpQGwm19/ld5/3xzXrSv94x+WlgMAnmDZreAaN26slStX6vjx49q7d68Mw9Cll16q2rVrW1USADuaPl1yOs3xgw9KISHW1gMAHmBZAMxXu3ZtXXPNNVaXAcCOTpyQ3nrLHIeESCNGWFkNAHiMJSeBAECVMHeuee9fSRo4UIqIsLYeAPAQAiAAe8rOlmbMKJx+4gnragEADyMAArCnd9+VUlLM8Z13Ss2bW1sPAHgQARCA/RiGNHVq4TTH/gGwGQIgAPtZvVravt0cX3ONees3ALARAiAA+3n55cLx6NGSw2FdLQBgAQIgAHvZtUv6+GNz3KiRdMcd1tYDABYgAAKwl7OP/Xv8cSnA8suhAoDHEQAB2EdqqrRwoTkODzfv/AEANkQABGAfM2ea1/+TpCFDpBo1LC0HAKxCAARgD6dOSa+9Zo4DA6Xhw62tBwAsRAAEYA8LFpghUDJv+xYTY2k5AGAlAiAA35eXJ02fXjj92GPW1QIAVYBPBcBZs2apSZMmCgkJUceOHbV58+bzWm/x4sVyOBzq1atX5RYIwBoffSTt22eO4+Kk1q2trQcALOYzAXDJkiWKj49XQkKCtm7dqtatW6t79+46cuRIqevt379fTz75pDp16uShSgF4lGFIiYmF0yNHWlcLAFQRPhMAp0+friFDhmjw4MFq1aqV5syZo9DQUM2fP9/tOnl5eerfv78mTpyoiy++2IPVAvCYTz+Vvv7aHLdpI91yi6XlAEBV4BNXQM3JydGWLVs0duzYgnl+fn6Ki4vTxo0b3a733HPPqV69err//vv1xRdfnPN9srOzlZ1/CQlJ6enpkiSn0ymn01kwPvsnCtGb0tEf9y6kN/4vvVTwl27uqFEycnMrsDLr8XvjHr1xz+69sevnPptPBMC0tDTl5eUpMjKyyPzIyEjt3r27xHXWr1+vt956S9u2bTvv90lMTNTEiROLzV+9erVCQ0OLzEtKSjrv17UbelM6+uNeWXsTvm+fuvyxzul69fRpSIi0cmVllGY5fm/cozfu2bU3mZmZVpdgOZ8IgGWVkZGhAQMGaN68eYqIiDjv9caOHav4+PiC6fT0dMXExKhbt24KDw+XZP5VkZSUpK5duyowMLDCa/dm9KZ09Me98vbG/777CsYhTz+tW3v2rIzyLMXvjXv0xj279yZ/D56d+UQAjIiIkL+/v1JTU4vMT01NVVRUVLHlf/75Z+3fv189z/qfgcvlkiQFBARoz549uuSSS4qtFxwcrODg4GLzAwMDi/0DKmkeTPSmdPTHvTL15sABackSc1yrlvwfeED+PtxXfm/cozfu2bU3dvzMf+YTJ4EEBQWpXbt2WrNmTcE8l8ulNWvWKDY2ttjyLVu21Pfff69t27YVPG6//XZ16dJF27ZtUwwXiAW83+TJUv5xPkOHStWrW1sPAFQhPrEFUJLi4+M1aNAgtW/fXh06dNDMmTN1+vRpDR48WJI0cOBANWjQQImJiQoJCdEVV1xRZP1atWpJUrH5ALzQgQPSW2+Z47Aw6cknra0HAKoYnwmAffr00dGjRzV+/HilpKSoTZs2WrVqVcGJIQcOHJCfn09s8ARwLi+9JOWf7Tt8uFSnjrX1AEAV4zMBUJKGDRumYcOGlfhccnJyqesuWLCg4gsC4HlHjkjz5pnj6tWlESMsLQcAqiI2iQHwLbNnS/nX63zwQakMZ/oDgF0QAAH4jpMnpZkzzbG/v/T445aWAwBVFQEQgO+YO9cMgZI0cKDUpIml5QBAVUUABOAbcnKk6dPNscMhjRplbT0AUIURAAH4hkWLpJQUc9yrl3TZZZaWAwBVGQEQgPfLy5OmTCmcPuuWjQCA4giAALzf0qXSnj3m+IYbpOuvt7YeAKjiCIAAvFtenjRhQuH0s89aVgoAeAsCIADv9u670o8/muPrr5duusnaegDACxAAAXivvDxp/PjC6WefNc8ABgCUigAIwHstWSL99JM5vvFGqWtXa+sBAC9BAATgnVwuadKkwulx49j6BwDniQAIwDstWiTt3m2OOfYPAMqEAAjA+zidRc/2TUhg6x8AlAEBEID3Wb5c2r/fHHfrJsXFWVkNAHgdAiAA72IY0ksvFU4/+aR1tQCAlyIAAvAuq1ZJX39tjq+6iq1/AFAOBEAA3sMwpMmTC6c59g8AyoUACMBrOD7+WFq/3pxo2VLq1cvSegDAWxEAAXgFR16e/EePLpwxYYLkx1cYAJQH354AvELDzz+XY88ecyI2Vurd29qCAMCLEQABVH1nzqjlokWF088/z7F/AHABCIAAqjy/t95S6NGj5sTNN0t//au1BQGAlyMAAqjajh+X3/PPF06fPQYAlAsBEEDV9vzzcvz+uyTJde+9Urt2FhcEAN6PAAig6tq3T3rlFUlSXmCg8iZNsrggAPANBEAAVdf48VJuriRpb69eUuPG1tYDAD6CAAigatqwQXr7bUmSUaeOGQABABWCAAig6nG5pEceKZx8+mnlVq9uYUEA4FsIgACqnv/8R9q+3Ry3bSvXWWEQAHDhCIAAqpasLGncuMLpKVOkgADr6gEAH0QABFC1TJsmHTxojm+7Tera1dp6AMAHEQABVB2pqdKLL5pjPz9z6x8AoMIRAAFUHaNHSxkZ5njIEOnyy62tBwB8FAEQQNWwe7f03/+a41q1pIkTLS0HAHwZARBA1TB6tHn5F0kaNUqKjLS2HgDwYQRAANZbt0766CNzHB0tDR9ubT0A4OMIgACs5XJJ8fGF05MmSVz0GQAqFQEQgLX+8x9p61Zz3Lq1NGiQtfUAgA0QAAFY59QpaezYwumXXpL8/a2rBwBsggAIwDrTpkkpKea4Vy8u+gwAHkIABGCNffukxERz7O/PRZ8BwIMIgAA8zzCkxx6TsrPN6eHDpRYtrK0JAGyEAAjA85KSpBUrzHH9+lJCgrX1AIDNEAABeFZOjjRyZOH0tGlSeLh19QCADREAAXjWyy9LO3ea42uukfr0sbYeALAhAiAAz9m3r/Aevw6H9Prrkh9fQwDgaXzzAvCckSOl06fN8dChUrt21tYDADZFAATgGevXSx9+aI6jogovAQMA8DgCIIDK53KZl3rJ9/zznPgBABYiAAKofDNnFt7v94oruN8vAFiMAAigcu3dKz39dOH0q69KAQHW1QMAIAACqESGIT3ySOEdP0aMkDp3trIiAIAIgAAq07//bd71Q5IaNpQmTbK2HgCAJAIggMpy5IgUH184PWuWFBZmXT0AgAIEQACV45FHpOPHzXG/ftLtt1tbDwCggE8FwFmzZqlJkyYKCQlRx44dtXnzZrfLzps3T506dVLt2rVVu3ZtxcXFlbo8gDJYvFh67z1zXKeOeb9fAECV4TMBcMmSJYqPj1dCQoK2bt2q1q1bq3v37jpy5EiJyycnJ6tv375au3atNm7cqJiYGHXr1k0HDx70cOWAjzl40LzLR77XXzcv/AwAqDJ8JgBOnz5dQ4YM0eDBg9WqVSvNmTNHoaGhmj9/fonLL1y4UI888ojatGmjli1b6s0335TL5dKaNWs8XDngQ1wu6e9/l06cMKd79zYfAIAqxScuxpWTk6MtW7Zo7NixBfP8/PwUFxenjRs3ntdrZGZmyul0qk6dOm6Xyc7OVnb+5SwkpaenS5KcTqecTmfB+OyfKERvSucL/fGbMUP+ycmSJKNhQ+W+/LJUAZ/HF3pTWeiNe/TGPbv3xq6f+2wOwzAMq4u4UIcOHVKDBg20YcMGxcbGFsx/6qmntG7dOm3atOmcr/HII4/ok08+0Y4dOxQSElLiMhMmTNDEiROLzV+0aJFCQ0PL/wEAH1D7xx91/dix8svLkyR9OWmS0q680uKqAKC4zMxM9evXTydPnlS4TW9L6RNbAC/UlClTtHjxYiUnJ7sNf5I0duxYxZ91WYv09PSCYwfzf4GcTqeSkpLUtWtXBQYGVnrt3oTelM6r+3P6tAKeekqOP8JfXny8OoweXWEv79W9qWT0xj16457de5O/B8/OfCIARkREyN/fX6mpqUXmp6amKuocB59PnTpVU6ZM0aeffqqrrrqq1GWDg4MVHBxcbH5gYGCxf0AlzYOJ3pTOK/szbpz044/muH17+U+ZIv9K+Axe2RsPoTfu0Rv37NobO37mP/OJk0CCgoLUrl27Iidw5J/QcfYu4T978cUXNWnSJK1atUrt27f3RKmA71m1yrzIsySFhEhvvy3x5QoAVZpPbAGUpPj4eA0aNEjt27dXhw4dNHPmTJ0+fVqDBw+WJA0cOFANGjRQYmKiJOmFF17Q+PHjtWjRIjVp0kQpKSmSpLCwMIVxtwLg/KSlSffdVzg9bZrUooV19QAAzovPBMA+ffro6NGjGj9+vFJSUtSmTRutWrVKkZGRkqQDBw7Iz69wg+frr7+unJwc3X333UVeJyEhQRMmTPBk6YD3Gj5cOnzYHHftWvT6fwCAKstnAqAkDRs2TMOGDSvxueQ/Lk2Rb//+/ZVfEODLli2TFi0yx7VqSQsWSA6HlRUBAM6TTxwDCMDDfvtNevDBwulXXpGio62rBwBQJgRAAGXjdJp39zh+3Jy+6y7z7h8AAK9BAARQNhMmSPl32GncWJo7l12/AOBlCIAAzt+aNdIfZ9IrIEBavFgq5faJAICqiQAI4Pz8+qvUr5+Uf/fISZOka6+1tiYAQLkQAAGcW2amdM890pEj5nTXrtJTT1lbEwCg3AiAAM7tscekTZvMcePG0jvvSH58fQCAt+IbHEDp/v1vaf58c1y9uvT++9JFF1lbEwDgghAAAbi3fXvR6/299prUtq119QAAKgQBEEDJ0tOlO++UcnLM6Ycekv7xD0tLAgBUDAIggOJcLmnQIOmXX8zpdu2kl1+2tiYAQIUhAAIo7umnpeXLzXHNmtLSpVJwsKUlAQAqDgEQQFEzZkgvvGCOHQ5p0SLp4outrQkAUKEIgAAKvfWWFB9fOP3qq9Ktt1pXDwCgUhAAAZhWrTJP9Mj3zDPSo49aVw8AoNIQAAGYF3m+6y4pL8+cHjHCvNUbAMAnEQABu9uxQ+rRw7zdmyT97W/StGnm8X8AAJ9EAATs7NAh876+v/9uTt94o3nSB7d5AwCfxrc8YFcpKdLNN0uHD5vT7dpJH34ohYRYWxcAoNIFWF0AAAukpkpduki7d5vTMTHSxx+b1/wDAPg8tgACdnPihNS9e9Hwl5Qk1a1raVkAAM9hCyBgJxkZ5m7f7dvN6ehoae1a6ZJLrK0LAOBRbAEE7OL0aem228xLvkjSRReZW/4IfwBgOwRAwA6OHZO6dZM+/9ycrl1b+uwzqVUra+sCAFiCXcCAr/v9d/OYvy1bzOnwcGn1aumqq6ytCwBgGbYAAr7sxAnpr38tDH/16pnH/LVvb2lZAABrEQABX3X0qLnb97vvzOmoKGnNGqltW2vrAgBYjl3AgC/au1e64w5p505zum5dad06qXlza+sCAFQJbAEEfM3KleZdPfLDX1SUebYv4Q8A8AcCIOArDEN64QXzUi/p6ea8li2lL76QWre2tjYAQJVCAAR8QW6uNGqUNGaMGQQlqVcv6auvpGbNLC0NAFD1cAwg4O0OHpTuuUfauLFw3nPPSePGSQ6HdXUBAKosAiDgzdavl/r0kQ4dMqcDAqRXXpGGDrW2LgBAlcYuYMAbGYYZ9Dp3Lgx/jRubd/og/AEAzoEtgIC3ycyUhgyRFi0qnHfDDdLixVL9+tbVBQDwGmwBBLzJ9u3SddcVDX+jR5sXeCb8AQDOE1sAAW9gGNLrr0sjR0o5Oea8sDBp/nzzBBAAAMqAAAhUdUeOSA8/LH3wQeG8yy+X3nlHuvJK6+oCAHgtdgEDVZVhSAsWSJddVjT8DR8ubdlC+AMAlBtbAIGqaP9+6YEHzGP78l10kfTGG9Jdd1lWFgDAN7AFEKhC/LOy5Ddlirl17+zw16eP9N13hD8AQIVgCyBQFeTmyvGvfylu9Gj5nzhROD8mRpozR7r1VstKAwD4HrYAAlZyuaQPP5TatFHAQw8pJD/8+flJDz0k/fAD4Q8AUOHYAghYwTCk5GTpqaekb74p8pTr9tvN3cCXXWZNbQAAn0cABDwpL09auVJ68UXzPr5ncXXooA133KGOo0bJLzDQogIBAHbALmDAE06elGbPllq2lG6/vWj4u/JK6cMPlffFF/r98sutqxEAYBtsAQQq065d0vTp5kWbT58u+lzLltLYsVL//pK/v+R0WlMjAMB2CIBARTt9Wlq2TPrPf6TPPiv+/I03Sk88IfXoYZ7sAQCAhxEAgYqQmyt9+qm5pe+994pv7atRQ/r736UHH5TatLGkRAAA8hEAgfLKzZXWrZPefde8VduRI8WXadZMGjJEeuQRKSzM8zUCAFACAiBQFikp0scfS6tXm2fzpqcXX6ZmTemee6R//EO67jrJ4fB4mQAAlIYACJTm1Cnpiy/MY/k+/VTatq3k5UJCzAs2//3v0i23mNMAAFRRBEDgbE6n9PXXUlKStHat9NVXUnZ2ycvWqiV17y797W/mCR3s4gUAeAkCIOzL6ZR+/FHauNG8G8e2bdL27dKZM+7XadfODH233ip16CBxwWYAgBciAMIeTp2Sfv5Z2rFD2rrV3LL39ddSTk7p6zVtKnXtaj5uuEGqV88z9QIAUIkIgPAdp0+bIW/vXmnfPmn//sLQd+DA+b3GJZdIHTtKN90kxcVJjRpVaskAAFjBpwLgrFmz9NJLLyklJUWtW7fWq6++qg4dOrhd/t1339Wzzz6r/fv369JLL9ULL7ygW2+91YMV47xkZkqHDpmPw4fNy62kpUlHj5qXYLn4Yum338znDOP8X/eSS8xdum3bmmfrtm4thYdX3ucAAKCK8JkAuGTJEsXHx2vOnDnq2LGjZs6cqe7du2vPnj2qV8Juuw0bNqhv375KTEzUbbfdpkWLFqlXr17aunWrrrjiCgs+gQ8zDDPEnTih0NRU6fvvpdRUc/fr8eNmmDtyxHwcPy6dOCEdO1Y4/vNFlf8sLa3052vUMO+326yZdNll5viaa9idCwCwLZ8JgNOnT9eQIUM0ePBgSdKcOXO0YsUKzZ8/X2PGjCm2/Msvv6ybb75Zo0aNkiRNmjRJSUlJeu211zRnzhyP1u6WYZiPvDzJ5TIfpY1zc4vOL+2n02k+cnMLf2Znm1vVatQww1lOjrnVLS/PvKzJmTNSVpb5M//57dul6tXNLWf5z58+bf7MyjKPvTt1SpIUKKlrZfarQwfz/rqXXGKGvcaNzWP46tfnWnwAAJzFJwJgTk6OtmzZorFjxxbM8/PzU1xcnDZu3FjiOhs3blR8fHyRed27d9fy5csrs9Tz06qVtHt32XZn+ppq1aTatc1LrdSpI0VHSw0bmmEuMlK66CJzC15UlDnN2bgAAJw3nwiAaWlpysvLU2RkZJH5kZGR2r17d4nrpKSklLh8SkqK2/fJzs5W9lnXhEv/4y4QTqdTTqezYHz2z/IIyMuTw0vDn+FwSKGhZoCrVk2qXl1GWJgUECDHjz/qYKtWirrkEjkuukiO48dlXHGFjMhIqW5dGfXqmcGuZk1z3bK4gH5XFRXxu+Or6I179MY9euOe3Xtj1899Np8IgJ6SmJioiRMnFpu/evVqhYaGFpmXlJRU7vdpHxGhUMOQ4ednBiqHwxz7+Zljh0OGv3+x5wx//6LL/bF+wbyzpl0BAeby/v5y/fHT8PeXIy9P/tnZyqpbV66AALkCAuTndMpZo4ZyQ0KUFxRkzs//GRgow99fucHBygsKkhEQUL7drVlZ0v/+Zz5s7kJ+d3wdvXGP3rhHb9yza28yMzOtLsFyPhEAIyIi5O/vr9TU1CLzU1NTFRUVVeI6UVFRZVpeksaOHVtkt3F6erpiYmLUrVs3hf9x9qjT6VRSUpK6du2qwPLulvTRM5ErpDc+jP64R2/cozfu0Rv37N6b9JLu424zPhEAg4KC1K5dO61Zs0a9evWSJLlcLq1Zs0bDhg0rcZ3Y2FitWbNGI0aMKJiXlJSk2NhYt+8THBys4ODgYvMDAwOL/QMqaR5M9KZ09Mc9euMevXGP3rhn197Y8TP/mU8EQEmKj4/XoEGD1L59e3Xo0EEzZ87U6dOnC84KHjhwoBo0aKDExERJ0vDhw3XjjTdq2rRp6tGjhxYvXqxvvvlGc+fOtfJjAAAAVDqfCYB9+vTR0aNHNX78eKWkpKhNmzZatWpVwYkeBw4ckJ+fX8Hy1113nRYtWqRx48bp6aef1qWXXqrly5dzDUAAAODzfCYAStKwYcPc7vJNTk4uNu+ee+7RPffcU8lVAQAAVC1+514EAAAAvoQACAAAYDMEQAAAAJshAAIAANgMARAAAMBmCIAAAAA2QwAEAACwGQIgAACAzRAAAQAAbMan7gTiaYZhSJLS09ML5jmdTmVmZio9PZ2bTf8JvSkd/XGP3rhHb9yjN+7ZvTf5/9/O//+4HREAL0BGRoYkKSYmxuJKAABAWWVkZKhmzZpWl2EJh2Hn+HuBXC6XDh06pBo1asjhcEgy/6qIiYnRr7/+qvDwcIsrrFroTenoj3v0xj164x69cc/uvTEMQxkZGYqOjpafnz2PhmML4AXw8/NTw4YNS3wuPDzclv+ozge9KR39cY/euEdv3KM37tm5N3bd8pfPnrEXAADAxgiAAAAANkMArGDBwcFKSEhQcHCw1aVUOfSmdPTHPXrjHr1xj964R2/ASSAAAAA2wxZAAAAAmyEAAgAA2AwBEAAAwGYIgAAAADZDAKwAx44dU//+/RUeHq5atWrp/vvv16lTp0pd/rHHHlOLFi1UrVo1NWrUSI8//rhOnjzpwao9o6y9kaS5c+eqc+fOCg8Pl8Ph0IkTJzxTrAfMmjVLTZo0UUhIiDp27KjNmzeXuvy7776rli1bKiQkRFdeeaVWrlzpoUo9ryy92bFjh+666y41adJEDodDM2fO9FyhFihLb+bNm6dOnTqpdu3aql27tuLi4s75e+bNytKb999/X+3bt1etWrVUvXp1tWnTRv/97389WK1nlfX7Jt/ixYvlcDjUq1evyi0QliIAVoD+/ftrx44dSkpK0v/7f/9Pn3/+uR588EG3yx86dEiHDh3S1KlT9cMPP2jBggVatWqV7r//fg9W7Rll7Y0kZWZm6uabb9bTTz/toSo9Y8mSJYqPj1dCQoK2bt2q1q1bq3v37jpy5EiJy2/YsEF9+/bV/fffr2+//Va9evVSr1699MMPP3i48spX1t5kZmbq4osv1pQpUxQVFeXhaj2rrL1JTk5W3759tXbtWm3cuFExMTHq1q2bDh486OHKK19Ze1OnTh0988wz2rhxo7777jsNHjxYgwcP1ieffOLhyitfWXuTb//+/XryySfVqVMnD1UKyxi4IDt37jQkGV9//XXBvI8//thwOBzGwYMHz/t1li5dagQFBRlOp7MyyrTEhfZm7dq1hiTj+PHjlVil53To0MF49NFHC6bz8vKM6OhoIzExscTle/fubfTo0aPIvI4dOxoPPfRQpdZphbL25myNGzc2ZsyYUYnVWetCemMYhpGbm2vUqFHD+Pe//11ZJVrmQntjGIZx9dVXG+PGjauM8ixVnt7k5uYa1113nfHmm28agwYNMu644w4PVAqrsAXwAm3cuFG1atVS+/btC+bFxcXJz89PmzZtOu/XOXnypMLDwxUQ4Du3Z66o3viCnJwcbdmyRXFxcQXz/Pz8FBcXp40bN5a4zsaNG4ssL0ndu3d3u7y3Kk9v7KIiepOZmSmn06k6depUVpmWuNDeGIahNWvWaM+ePbrhhhsqs1SPK29vnnvuOdWrV88n90ahON9JGxZJSUlRvXr1iswLCAhQnTp1lJKScl6vkZaWpkmTJp1z16i3qYje+Iq0tDTl5eUpMjKyyPzIyEjt3r27xHVSUlJKXN7Xelee3thFRfRm9OjRio6OLvbHhLcrb29OnjypBg0aKDs7W/7+/po9e7a6du1a2eV6VHl6s379er311lvatm2bBypEVcAWQDfGjBkjh8NR6qMi/ueUnp6uHj16qFWrVpowYcKFF+4BnuoNgAszZcoULV68WB988IFCQkKsLqdKqFGjhrZt26avv/5a//znPxUfH6/k5GSry7JURkaGBgwYoHnz5ikiIsLqcuAhbAF044knntA//vGPUpe5+OKLFRUVVeyg2tzcXB07duycB6dnZGTo5ptvVo0aNfTBBx8oMDDwQsv2CE/0xtdERETI399fqampReanpqa67UVUVFSZlvdW5emNXVxIb6ZOnaopU6bo008/1VVXXVWZZVqivL3x8/NTs2bNJElt2rTRrl27lJiYqM6dO1dmuR5V1t78/PPP2r9/v3r27Fkwz+VySTL32uzZs0eXXHJJ5RYNj2MLoBt169ZVy5YtS30EBQUpNjZWJ06c0JYtWwrW/eyzz+RyudSxY0e3r5+enq5u3bopKChIH330kVf9dV7ZvfFFQUFBateundasWVMwz+Vyac2aNYqNjS1xndjY2CLLS1JSUpLb5b1VeXpjF+XtzYsvvqhJkyZp1apVRY7B9SUV9XvjcrmUnZ1dGSVapqy9admypb7//ntt27at4HH77berS5cu2rZtm2JiYjxZPjzF6rNQfMHNN99sXH311camTZuM9evXG5deeqnRt2/fgud/++03o0WLFsamTZsMwzCMkydPGh07djSuvPJKY+/evcbhw4cLHrm5uVZ9jEpR1t4YhmEcPnzY+Pbbb4158+YZkozPP//c+Pbbb43ff//dio9QYRYvXmwEBwcbCxYsMHbu3Gk8+OCDRq1atYyUlBTDMAxjwIABxpgxYwqW//LLL42AgABj6tSpxq5du4yEhAQjMDDQ+P777636CJWmrL3Jzs42vv32W+Pbb7816tevbzz55JPGt99+a/z0009WfYRKU9beTJkyxQgKCjKWLVtW5LslIyPDqo9Qacram8mTJxurV682fv75Z2Pnzp3G1KlTjYCAAGPevHlWfYRKU9be/BlnAfs+AmAF+P33342+ffsaYWFhRnh4uDF48OAiX7b79u0zJBlr1641DKPw8iYlPfbt22fNh6gkZe2NYRhGQkJCib3517/+5fkPUMFeffVVo1GjRkZQUJDRoUMH46uvvip47sYbbzQGDRpUZPmlS5cazZs3N4KCgozLL7/cWLFihYcr9pyy9Cb/9+bPjxtvvNHzhXtAWXrTuHHjEnuTkJDg+cI9oCy9eeaZZ4xmzZoZISEhRu3atY3Y2Fhj8eLFFlTtGWX9vjkbAdD3OQzDMDy3vREAAABW4xhAAAAAmyEAAgAA2AwBEAAAwGYIgAAAADZDAAQAALAZAiAAAIDNEAABAABshgAIAABgMwRAAChFcnKy2rZtq+DgYDVr1kwLFiywuiQAuGAEQABwY9++ferRo4e6dOmibdu2acSIEXrggQf0ySefWF0aAFwQbgUHwLbmzp2rCRMm6LfffpOfX+Hfw3fccYcuuugi1a1bVytWrNAPP/xQ8Ny9996rEydOaNWqVVaUDAAVgi2AAGzrnnvu0e+//661a9cWzDt27JhWrVql/v37a+PGjYqLiyuyTvfu3bVx40ZPlwoAFYoACMC2ateurVtuuUWLFi0qmLds2TJFRESoS5cuSklJUWRkZJF1IiMjlZ6erqysLE+XCwAVhgAIwNb69++v9957T9nZ2ZKkhQsX6t577y2ySxgAfA3fcABsrWfPnjIMQytWrNCvv/6qL774Qv3795ckRUVFKTU1tcjyqampCg8PV7Vq1awoFwAqRIDVBQCAlUJCQnTnnXdq4cKF2rt3r1q0aKG2bdtKkmJjY7Vy5coiyyclJSk2NtaKUgGgwrAFEIDt9e/fXytWrND8+fMLtv5J0sMPP6xffvlFTz31lHbv3q3Zs2dr6dKlGjlypIXVAsCF4zIwAGzP5XKpYcOGOnz4sH7++WddfPHFBc8lJydr5MiR2rlzpxo2bKhnn31W//jHP6wrFgAqAAEQAADAZtgFDAAAYDMEQAAAAJshAAIAANgMARAAAMBmCIAAAAA2QwAEAACwGQIgAACAzRAAAQAAbIYACAAAYDMEQAAAAJshAAIAANgMARAAAMBm/j9b+zGKP0W+IgAAAABJRU5ErkJggg==",
      "text/plain": [
       "class=Graph name=v0 CDF implementation=class=GraphImplementation name=v0 CDF title=v0 CDF xTitle=v0 yTitle=CDF axes=ON grid=ON legendposition=topleft legendFontSize=10 drawables=[class=Drawable name=v0 CDF implementation=class=Staircase name=v0 CDF pattern=s derived from class=DrawableImplementation name=v0 CDF legend=v0 CDF data=class=Sample name=Unnamed implementation=class=SampleImplementation name=Unnamed size=99968 dimension=2 data=[[-0.175795,0.00031],[-0.175113,0.00032],[-0.174854,0.00033],...,[0.455322,0.99996],[0.456101,0.99997],[0.458048,0.99997]] color=red isColorExplicitlySet=true fillStyle=solid lineStyle=solid pointStyle=plus lineWidth=2]"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Visualize the clearance CDF\n",
    "clearance_distribution.drawCDF()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f822b3a8-86e0-4d50-9278-cca3df39e154",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Raw Cell Format",
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}