{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "8ef12ff3",
   "metadata": {},
   "source": [
    "# Fragility curves\n",
    "\n",
    "### Failure probability calculations with a fragility curve\n",
    "\n",
    "In this example, we will demonstrate how to calculate the failure probability of a levee using a fragility curve.\n",
    "\n",
    "Fragility curve is a function that describes the relation between the load imposed on the levee and the corresponding  (conditional) failure probability. Typically a fragility curve is expressed by the following relation: $h \\rightarrow P(Z<0|h)$. By integrating the fragility curve with the load statistics, the failure probability can be estimated. The goal is to derive the following probability: \n",
    "\n",
    "$P(Z<0) = \\int P(Z<0 | h)\\cdot f(h) dh$\n",
    "\n",
    "First, let's import the necessary packages:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "4989b056",
   "metadata": {},
   "outputs": [],
   "source": [
    "from probabilistic_library import ReliabilityProject, DistributionType, ReliabilityMethod, FragilityCurve, FragilityValue, Stochast\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a853414c",
   "metadata": {},
   "source": [
    "We consider the Hunt's limit state function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "16a403f5",
   "metadata": {},
   "outputs": [],
   "source": [
    "from utils.models import hunt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "456f8cb5",
   "metadata": {},
   "source": [
    "We define a reliability project using the `ReliabilityProject()` class and refer to the Hunt's limit state function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "03ef4c74",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Model hunt:\n",
      "Input parameters:\n",
      "  t_p\n",
      "  tan_alpha\n",
      "  h_s\n",
      "  h_crest\n",
      "  h\n",
      "Output parameters:\n",
      "  Z\n"
     ]
    }
   ],
   "source": [
    "project = ReliabilityProject()\n",
    "project.model = hunt\n",
    "\n",
    "project.model.print()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ecde9524",
   "metadata": {},
   "source": [
    "We define the associated variables, except for the water level variable, $h$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "35b70f41",
   "metadata": {},
   "outputs": [],
   "source": [
    "project.variables[\"t_p\"].distribution = DistributionType.log_normal\n",
    "project.variables[\"t_p\"].mean = 6\n",
    "project.variables[\"t_p\"].deviation = 2\n",
    "\n",
    "project.variables[\"tan_alpha\"].distribution = DistributionType.deterministic\n",
    "project.variables[\"tan_alpha\"].mean = 0.333333\n",
    "\n",
    "project.variables[\"h_s\"].distribution = DistributionType.log_normal\n",
    "project.variables[\"h_s\"].mean = 3\n",
    "project.variables[\"h_s\"].deviation = 1\n",
    "\n",
    "project.variables[\"h_crest\"].distribution = DistributionType.log_normal\n",
    "project.variables[\"h_crest\"].mean = 10\n",
    "project.variables[\"h_crest\"].deviation = 0.05"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f56ba43c",
   "metadata": {},
   "source": [
    "We choose the `form` calculation method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "5f8e71c8",
   "metadata": {},
   "outputs": [],
   "source": [
    "project.settings.reliability_method = ReliabilityMethod.form\n",
    "project.settings.relaxation_factor = 0.15\n",
    "project.settings.maximum_iterations = 50\n",
    "project.settings.epsilon_beta = 0.01"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7ce6f4b5",
   "metadata": {},
   "source": [
    "Next, we construct the fragility curve using the `FragilityCurve()` class, providing the water level variable $h$. Reliability calculations are performed for each value of the water level variable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "7a390139",
   "metadata": {},
   "outputs": [],
   "source": [
    "fragility_curve = FragilityCurve()\n",
    "fragility_curve.name = \"conditional\"\n",
    "\n",
    "fc_pf = []\n",
    "h = np.arange(-2.0, 5.0, 0.25)\n",
    "for val in h:\n",
    "\n",
    "        project.variables[\"h\"].distribution = DistributionType.deterministic\n",
    "        project.variables[\"h\"].mean = val\n",
    "        project.run()\n",
    "        dp = project.design_point\n",
    "        \n",
    "        value = FragilityValue()\n",
    "        value.x = val\n",
    "        value.reliability_index = dp.reliability_index\n",
    "        value.design_point = dp\n",
    "        \n",
    "        fragility_curve.fragility_values.append(value)\n",
    "        fc_pf.append(dp.probability_failure)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1cd49c68",
   "metadata": {},
   "source": [
    "The calculations above result in a fragility curve:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "afa76689",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": [
     "gallery",
     "reliability"
    ]
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGwCAYAAABB4NqyAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABLB0lEQVR4nO3de1zUVf4/8NfMMMxwHUUuA3gBFTU08ArhJS1RtDIt19QyL19Xv7na6o/MsjWRale3rOxiWpamXzWr3bK1C0mUVhuCSmqI91BU7hIXQWCY+fz+QGZEZmCAYT5zeT0fDx8rn/nM4T0nWl6dcz7nSARBEEBERETkRKRiF0BERERkbQxARERE5HQYgIiIiMjpMAARERGR02EAIiIiIqfDAEREREROhwGIiIiInI6L2AXYKp1Oh9zcXHh5eUEikYhdDhEREZlBEARUVFQgKCgIUqnpcR4GIBNyc3PRrVs3scsgIiKiNrh8+TK6du1q8nUGIBO8vLwA1Hegt7e3xdrVaDTYv38/xo8fD7lcbrF27RH7woB9UY/9YMC+MGBfGLAvDEz1RXl5Obp166b/PW4KA5AJDdNe3t7eFg9A7u7u8Pb25g8v+0KPfVGP/WDAvjBgXxiwLwxa6ouWlq9wETQRERE5HQYgIiIicjoMQEREROR0GICIiIjI6TAAERERkdNhACIiIiKnwwBERERETocBiIiIiJwOAxARERE5He4ETURERFaj1QlIzy5BYUU1/L2UiAr1gUxq/UPHGYCIiIjIKpIy85C4Lwt5ZdX6a4EqJRImhWPCgECr1sIpMCIiIupwSZl5WLQzo1H4AYD8smos2pmBpMw8q9bDAEREREQdSqsTkLgvC4KR1xquJe7LglZn7I6OwQBEREREHSo9u6TJyM+tBAB5ZdVIzy6xWk0MQERERNShCitMh5+23GcJDEBERETUofy9lBa9zxIYgIiIiKhDRYX6IFClhKmH3SWofxosKtTHajUxABEREVGHkkklSJgUDgBNQlDD1wmTwq26HxADEBEREXW4CQMCMeuuHvBUNN6CUK1SYtOswVbfB4gbIRIREZFVHLn0Bypq6rBodC/0C/TiTtBERETk2K78UYVTeeWQSoAFd/eEj4erqPVwCoyIiIg63HdZBQCAoT18RA8/AAMQERERWcF3pwoBAOPCA0SupB4DEBEREXWo8moNDv1+DQAQywBEREREzuDgmSLU6QT08vNAqK+H2OUAYAAiIiKiDlZYUQN3V5nNjP4AfAqMiIiIOtj8kaF4LLo7ajQ6sUvRYwAiIiKiDqeUy6CUy8QuQ49TYERERNRhyqo0YpdgFAMQERERdQhBEPDA2z9h7KsHcK6gQuxyGuEUGBEREXWIswXXcbnkBlxdpAju7CZ2OY1wBIiIiIg6xHen6nd/HtnbF+6utjXmwgBEREREHSL55vEXsXfYzuPvDRiAiIiIyOIKK6px7HIpAGDsHf7iFmMEAxARERFZ3Pc3z/6K7NYJAd5KkatpigGIiIiILK5h/c84Gxz9AfgUGBEREXWAWXf1gK+nAnH91WKXYhQDEBEREVncmL7+GNPXNkd/AE6BERERkRNiACIiIiKL0eoEbPjuLH7N+QM6nSB2OSYxABEREZHFHLv8BzZ8dw6zt6ZDKzAAERERkRNIzqp//P2evv6Qy2w3ZthuZURERGR39I+/h9ve7s+3YgAiIiIii8gursT5wutwkUowuq+f2OU0iwGIiIiILCLl5ujPXT27wFspF7ma5jEAERERkUXs1x9+arv7/zRgACIiIqJ2q9Zo8XvRdQDAWBs8/f123AmaiIiI2k0pl+HQyrHIzC1HNx93sctpEUeAiIiIyCJcZFIM7NZJ7DLMwgBERERE7aLTCRBseNNDY2wiAG3cuBEhISFQKpWIjo5Genq6yXu3bNmCUaNGoXPnzujcuTNiY2Ob3C8IAlavXo3AwEC4ubkhNjYW586d6+iPQURE5JR+Ol+M4eu+x+vJZ8UuxWyiB6CPP/4Y8fHxSEhIQEZGBiIjIxEXF4fCwkKj9x84cAAzZ87EDz/8gNTUVHTr1g3jx4/H1atX9fe8/PLLePPNN7F582akpaXBw8MDcXFxqK6uttbHIiIichrJWfnIK6tG0fUasUsxm+gB6LXXXsOCBQswb948hIeHY/PmzXB3d8fWrVuN3r9r1y785S9/wcCBA9GvXz+8//770Ol0SElJAVA/+rNhwwasWrUKkydPRkREBHbs2IHc3Fzs3bvXip+MiIjI8QmCgO9uHn8xzg6e/mog6lNgtbW1OHr0KFauXKm/JpVKERsbi9TUVLPaqKqqgkajgY+PDwAgOzsb+fn5iI2N1d+jUqkQHR2N1NRUzJgxw2g7NTU1qKkxJNfy8nIAgEajgUajafVnM6WhLUu2aa/YFwbsi3rsBwP2hQH7wsAW++Jkbjnyy6vhJpdiWHdvq9Vmqi/M/f6iBqDi4mJotVoEBDROjAEBATh9+rRZbTzzzDMICgrSB578/Hx9G7e32fCaMWvXrkViYmKT6/v374e7u+Uf50tOTrZ4m/aKfWHAvqjHfjBgXxiwLwxsqS++uSwFIEWYVx1Skr+1+ve/vS+qqqrMep9d7wO0bt067NmzBwcOHIBSqWxXWytXrkR8fLz+6/Lycv36Im9v7/aWqqfRaJCcnIxx48ZBLrftbcI7GvvCgH1Rj/1gwL4wYF8Y2GJfbN6YCqACj425E/cNCrba9zXVFw0zOC0RNQD5+vpCJpOhoKCg0fWCggKo1epm37t+/XqsW7cO3333HSIiIvTXG95XUFCAwMDARm0OHDjQZHsKhQIKhaLJdblc3iE/ZB3Vrj1iXxiwL+qxHwzYFwbsCwNb6YurpTdwKr8CUgkQGx4oSk2394W5NYi6CNrV1RVDhgzRL2AGoF/QHBMTY/J9L7/8Ml588UUkJSVh6NChjV4LDQ2FWq1u1GZ5eTnS0tKabZOIiIhaRwJg3ogQTIoMQhfPpoMItkz0KbD4+HjMmTMHQ4cORVRUFDZs2IDKykrMmzcPADB79mwEBwdj7dq1AIB//vOfWL16NXbv3o2QkBD9uh5PT094enpCIpFg2bJleOmllxAWFobQ0FA8//zzCAoKwpQpU8T6mERERA4nqJMbEib1F7uMNhE9AE2fPh1FRUVYvXo18vPzMXDgQCQlJekXMefk5EAqNQxUbdq0CbW1tfjTn/7UqJ2EhASsWbMGALBixQpUVlZi4cKFKC0txciRI5GUlNTudUJERETkGEQPQACwZMkSLFmyxOhrBw4caPT1xYsXW2xPIpHghRdewAsvvGCB6oiIiOh2mVfLUHZDg6hQH8hlom8r2Gr2VzERERGJ7v2ffsdj76fZ1fEXt2IAIiIiolbRaHX4/nT97s/39PMXuZq2YQAiIiKiVjly8Q+UV9ehs7scg7t3FrucNrGJNUBERERk+7Q6AenZJXj3xwsAgHv6+kMmlYhcVdswABEREVGLkjLzkLgvC3ll1fpr358pRFJmHiYMCGzmnbaJU2BERETUrKTMPCzamdEo/ABAWZUGi3ZmICkzT6TK2o4BiIiIiEzS6gQk7suCYOS1hmuJ+7Kg1Rm7w3YxABEREZFJ6dklTUZ+biUAyCurRnp2ifWKsgAGICIiIjKpsMJ0+GnLfbaCAYiIiIhM8vcy7xgpc++zFQxAREREZFJUqA8CVUqYethdAiBQpURUqI81y2o3BiAiIiIySSaVIGFSuNFF0A2hKGFSuN3tB8QARERERM2aMCAQ9/T1a3JdrVJi06zBdrkPEDdCJCIiomZptDqcuFIGAHgmri+COrvB36t+2sveRn4aMAARERFRs34+X4xrlbXo4uGKP9/dE3KZ/U8g2f8nICIiog71xa9XAQAPRAQ6RPgBGICIiIioGVW1ddifVQAAmDwoWORqLIdTYERERGRSnU7AE6N74eilPzCoWyexy7EYBiAiIiIyyVspx1/HholdhsVxCoyIiIicDgMQERERGfXD6UJ8/VseqjVasUuxOAYgIiIiMuqNlHP4y64MfJSeI3YpFscARERERE1cLK7EsculkEqA+yPsb6fnljAAERERURNfHMsFAIwM87O7k97NwQBEREREjQiCgL3H6jc/nDIwSORqOgYDEBERETVy4koZsosroZRLMb6/WuxyOgQDEBERETXSMPozPlwNT4VjbhnIAERERESNXC6pAgBMGeSY018Ad4ImIiKi27w/Zxh+L7qObj7uYpfSYRiAiIiIqImefp5il9ChOAVGREREAIBqjRZlNzRil2EVDEBEREQEAPj2ZD6GvfQd1vznpNildDgGICIiIgJQv/lhrVYHL6Xjr5BhACIiIiJcu16Dg2eLAACTBwaLXE3HYwAiIiIifPVbHrQ6AXcGq9Db37EXQAMMQERERARg76/1mx9OdtCjL27HAEREROTkcq5VISOn/uT3ByMZgIiIiMgJfHHz6IsRvX3h7+14J78b4/jLvImIiKhZ04d1g7vCBT19PcQuxWoYgIiIiJycv7cS80eGil2GVXEKjIiIiJwOAxAREZGT0uoELN6VgU+OXEZNnVbscqyKAYiIiMhJ/XKhGF/9lod/fH0KEkjELseqGICIiIic1N5fcwEA998ZCFcX54oEzvVpiYiICABwo1aLpMw8AMCUQY5/9MXtGICIiIic0HenClBZq0XXzm4Y0r2z2OVYHQMQERGRE2rY/HDywCBIpc61/gdgACIiInI6JZW1OHCm/uT3KU5w8rsx3AiRiIjIyZRU1iK6pw/Kb9QhLMBL7HJEwQBERETkJLQ6AenZJSisqMaSe8IQ2U0ldkmiYQAiIiJyAkmZeUjcl4W8smr9tUCVEgmTwjFhQKCIlYmDa4CIiIgcXFJmHhbtzGgUfgAgv6wai3Zm6B+HdyYMQERERA5MqxOQuC8LgpHXGq4l7suCVmfsDsfFAEREROTA0rNLmoz83EoAkFdWjfTsEusVZQMYgIiIiBxYYYXp8NOW+xwFAxAREZED8/dSWvQ+R8EARERE5MCiQn0QqDIdbiSofxosKtTHekXZAAYgIiIiByaTSrD6gXCjrzUcgJEwKRwyJzsOgwGIiIjIwfVVe8FVJsXtEUetUmLTrMFOuQ8QN0IkIiJycD39PHF4VSx+u1wKmUyKwopq+HvVT3s528hPAwYgIiIiJ6Byk2NkHz+xy7AZnAIjIiJyYM72eLu5GICIiIgcVE2dFve98RMefPtnXC29IXY5NoVTYERERA4qKTMfxddrIZNK4O+lELscm2ITI0AbN25ESEgIlEoloqOjkZ6ebvLekydPYurUqQgJCYFEIsGGDRua3LNmzRpIJJJGf/r169eBn4CIiMj2bP/lIgDgsegekMts4le+zRC9Nz7++GPEx8cjISEBGRkZiIyMRFxcHAoLC43eX1VVhZ49e2LdunVQq9Um2+3fvz/y8vL0f37++eeO+ghEREQ2J/NqGTJySiGXSTAjqpvY5dgc0afAXnvtNSxYsADz5s0DAGzevBlfffUVtm7dimeffbbJ/cOGDcOwYcMAwOjrDVxcXJoNSLerqalBTU2N/uvy8nIAgEajgUajMbudljS0Zck27RX7woB9UY/9YMC+MGBfGLSmLz78bzYAIC48AJ2VMofrP1N9Ye7nFDUA1dbW4ujRo1i5cqX+mlQqRWxsLFJTU9vV9rlz5xAUFASlUomYmBisXbsW3bt3N3n/2rVrkZiY2OT6/v374e7u3q5ajElOTrZ4m/aKfWHAvqjHfjBgXxiwLwxa6otKDfDFrzIAEvTSXcHXX1+xTmEiuL0vqqqqzHqfqAGouLgYWq0WAQEBja4HBATg9OnTbW43OjoaH374Ifr27Yu8vDwkJiZi1KhRyMzMhJeXl9H3rFy5EvHx8fqvy8vL0a1bN4wfPx7e3t5truV2Go0GycnJGDduHORyucXatUfsCwP2RT32gwH7woB9YWBuX/zfoRxohNO4Q+2FxdPvgkTieJsdmuqLhhmclog+BdYRJk6cqP97REQEoqOj0aNHD3zyySeYP3++0fcoFAooFE1XyMvl8g75F66j2rVH7AsD9kU99oMB+8KAfWHQUl/MHh6KYB8PKFykcHV1tWJl1nd7X5j7MyJqAPL19YVMJkNBQUGj6wUFBa1av9OSTp06oU+fPjh//rzF2iQiIrJVLjIp4vpb7veoIxL1KTBXV1cMGTIEKSkp+ms6nQ4pKSmIiYmx2Pe5fv06Lly4gMBA5zvsjYiInItOJ4hdgl0Q/TH4+Ph4bNmyBdu3b8epU6ewaNEiVFZW6p8Kmz17dqNF0rW1tTh27BiOHTuG2tpaXL16FceOHWs0urN8+XIcPHgQFy9exC+//IKHHnoIMpkMM2fOtPrnIyIispaLxZUY8c/v8WbKOQgCg1BzRF8DNH36dBQVFWH16tXIz8/HwIEDkZSUpF8YnZOTA6nUkNNyc3MxaNAg/dfr16/H+vXrMXr0aBw4cAAAcOXKFcycORPXrl2Dn58fRo4ciUOHDsHPj4fAERGR49p56BLyyqqRkfOHQy58tiTRAxAALFmyBEuWLDH6WkOoaRASEtJiqt2zZ4+lSiMiIrILN2q1+OTIZQDA7JgeIldj+0SfAiMiIqL2++LYVZRX16G7jztG9/EXuxybxwBERERk5wRBwPbUSwCAWXd1h0zK6a+WMAARERHZuaOX/sCpvHIoXKR4ZCjP/TIHAxAREZGdaxj9mTwwCJ3cHXvjQ0uxiUXQRERE1HZPjO4JD1cZZt3Fxc/mYgAiIiKyc/2DVFg3NULsMuwKp8CIiIjI6TAAERER2an9J/Ox/NPjOJlbJnYpdodTYERERHbqg5+zkZZdArW3Ev2DVGKXY1c4AkRERGSHzuRXIC27BDKpBI9Gdxe7HLvDAERERGSHdqReBACMuyMAQZ3cxC3GDjEAERER2ZmKag0+//UqAGD2cD763hYMQERERHbm82N5qKrVIszfEzE9u4hdjl1iACIiIrIjggDsSssBUH/qu0TCc7/agk+BERER2QGtTkBadgkOF0kwqHsnSCUSPDS4q9hl2S0GICIiIhuXlJmHxH1ZyCurBiADkAu1SomfzxVhwoBAscuzS5wCIyIismFJmXlYtDPjZvgxKCirxqKdGUjKzBOpMvvGAERERGSjtDoBifuyIBh5reFa4r4saHXG7qDmMAARERHZqPTskiYjP7cSAOSVVSM9u8R6RTmIdq0BysnJwaVLl1BVVQU/Pz/0798fCoXCUrURERE5tcIK0+GnLfeRQasD0MWLF7Fp0ybs2bMHV65cgSAYht1cXV0xatQoLFy4EFOnToVUygEmIiKitvL3Ulr0PjJoVUL561//isjISGRnZ+Oll15CVlYWysrKUFtbi/z8fHz99dcYOXIkVq9ejYiICBw+fLij6iYiInJ4UaE+CFSZDjcSAIEqJaJCfaxXlINo1QiQh4cHfv/9d3Tp0nTXSX9/f9x777249957kZCQgKSkJFy+fBnDhg2zWLFERETORCaV4Ln7+uHJj441ea1h+8OESeGQSbkZYmu1KgCtXbvW7HsnTJjQ6mKIiIiosdIbdQAAqQS49WEvtUqJhEnh3AeojbgRIhERkQ3r6euBiK4qPDw4GL193bH/pzSMHxWNmN7+HPlphzYHoIKCAixfvhwpKSkoLCxstBgaALRabbuLIyIicnYjevvii8UjoBMAnbYO104JiA71YfhppzYHoLlz5yInJwfPP/88AgMDeRgbERFRB5FIJJBJAB3HFiymzQHo559/xk8//YSBAwdasBwiIiICgM9/vYKrf9zAnOEh8FLKxS7H4bQ5AHXr1q3JtBcRERG1X22dDuu/PYurpTfgpZRjzvAQsUtyOG3eqXDDhg149tlncfHiRQuWQ0RERP/OuIKrpTfg56XA9GHdxC7HIbVqBKhz586N1vpUVlaiV69ecHd3h1zeeHiupITnkhAREbVWbZ0Ob39/HgDwxOheUMplIlfkmFoVgDZs2NBBZRAREREAfHZz9MfXU4HHoruLXY7DalUAmjNnTkfVQURE5PQ0Wh3e/qFh9KcnR386EE8rJSIishGfZVzBlT8aRn96iF2OQ2vVCNC8efPatN/PlClT8OCDD7b6fURERM5kWIgPHhoUjAHBKri5cvSnI7UqAIWEhLTpm3Tq1KlN7yMiInImPf088fr0gWKX4RRaFYASEhI6qg4iIiIiq+EaICIiIpF9lnEFy/b8igtF18UuxWm0aSfo4uJibN26FampqcjPzwcAqNVqDB8+HHPnzoWfn59FiyQiInJUdVod3kg5h0vXqtAv0Bu9RnuKXZJTaPUI0OHDh9GnTx+8+eabUKlUuPvuu3H33XdDpVLhzTffRL9+/XDkyJGOqJWIiMjh7D2Wi0vXquDj4YrH7+KTX9bS6hGgJ598EtOmTcPmzZubPBEmCAKeeOIJPPnkk0hNTbVYkURERI6oTqvD29+fAwAsGNUTHoo2H9FJrdTqnj5+/Dg+/PBDo4/DSyQS/L//9/8waNAgixRHRETkyP5zPBcXr1Whs7scs2M4+mNNrZ4CU6vVSE9PN/l6eno6AgIC2lUUERGRo6vT6vDWzTO/FtzN0R9ra3VvL1++HAsXLsTRo0cxduxYfdgpKChASkoKtmzZgvXr11u8UCIiIkfy5Yk8ZBdXopO7HLNjQsQux+m0OgAtXrwYvr6+eP311/HOO+9Aq9UCAGQyGYYMGYIPP/wQjzzyiMULJSIiciRj7/DH8vF94KWUw5OjP1bXph6fPn06pk+fDo1Gg+LiYgCAr68v5HK5RYsjIiJyVF5KOZbcGyZ2GU6rXZFTLpcjMDDQUrUQERE5PJ1OgESCNp2tSZbTpp2gc3Nz8dRTT6G6urrR9bq6Ojz99NO4ePGiJWojIiJyOPtO5GLyxv/ix7NFYpfi1NoUgAIDA/HJJ5/g888/b3T9m2++wfbt29G1a1eLFEdEROQItDoBqReu4fNfr+Kf35zGiStlOH65VOyynFqbpsAkEglmzZqFHTt2YObMmfrr//d//4cZM2bAxYWLuYiIiAAgKTMPifuykFdmmDWRAAju7CZeUdT2w1DnzJmDlJQUFBQUAADKysrw5ZdfYu7cuZaqjYiIyK4lZeZh0c6MRuEHAAQAT31yHEmZeeIURm0PQP369cOQIUOwa9cuAMCnn36KXr16YfDgwRYrjoiIyF5pdQIS92VBaOaexH1Z0Oqau4M6SpsDEADMnj0bO3bsAFA//TVnzhyLFEVERGTv0rNLmoz83EoAkFdWjfTsEusVRXrtCkAzZ87EmTNn8MUXX+DQoUN4/PHHLVUXERGRXSusMB1+2nIfWVa7AlCnTp3w4IMPYv78+YiNjeUZYERERDf5eykteh9ZVrsCEFA/DVZSUsLpLyIioltEhfogUKWEqe0OJQACVUpEhfpYsyy6qd0B6L777kN2djamTp1qiXqIiIgcgkwqwf+O7gkATUJQw9cJk8Ihk3JHaDG0OwBJJBL06NEDMpnMEvUQERE5BK1OwCeHr8BD4QIfD9dGr6lVSmyaNRgTBvA4KbG0esfC+fPnY+HChYiOjjb6+h9//IGpU6fi+++/b3dxRERE9mp3eg6y8srhrXTBN0tH4UJRJQorquHvVT/txZEfcbV6BGjbtm0YM2YMtm3bZvT12tpaHDx4sN2FERER2auSylqs//YMAOCp8X3h761ETK8umDwwGDG9ujD82IA2TYGtWLEC//u//4ulS5dCp9NZuiYiIiK7tn7/GZTd0KCf2guPRXcXuxwyok0BaPHixUhOTsaePXswfvx4lJRwEyciIiIA+O1KGT5KzwEAvDB5AFxk7V5uSx2gzf9URo8ejfT0dFy7dg3Dhg1DZmamJesiIiKyOzqdgNX/yYQgAJMHBvERdxvWrljao0cP/PLLL4iKisLw4cPx2WeftamdjRs3IiQkBEqlEtHR0UhPTzd578mTJzF16lSEhIRAIpFgw4YN7W6TiIjIEmrqdAjz94SnwgXP3XeH2OVQM1odgCSSxgu33Nzc8NFHH+Fvf/sbZsyYgcTExFa19/HHHyM+Ph4JCQnIyMhAZGQk4uLiUFhYaPT+qqoq9OzZE+vWrYNarbZIm0RERJbg5irDy3+KxMGnxyDAmzs827JWByBBMH5q7TPPPIMvvvgCe/bsaVV7r732GhYsWIB58+YhPDwcmzdvhru7O7Zu3Wr0/mHDhuGVV17BjBkzoFAoLNImERGRJXXxNP77iWxHq/cB2rZtG1QqldHXJk6ciLS0NHz00UdmtVVbW4ujR49i5cqV+mtSqRSxsbFITU1tbWntarOmpgY1NTX6r8vLywEAGo0GGo2mTbUY09CWJdu0V+wLA/ZFPfaDAfvCwB764lzBdbyech7PTOiDHj7uHfZ97KEvrMVUX5jbN60OQC2d+RUWFobVq1eb1VZxcTG0Wm2TQ1QDAgJw+vTp1pbWrjbXrl1rdPpu//79cHe3/A9zcnKyxdu0V+wLA/ZFPfaDAfvCwFb7QhCAjVlSnCuXorgwH3P7dPz2MLbaF2K4vS+qqqrMel+rAtC6deuwdOlSuLm5tXhvWloaiouLcf/997fmW4hm5cqViI+P139dXl6Obt26Yfz48fD29rbY99FoNEhOTsa4ceMgl8st1q49Yl8YsC/qsR8M2BcGtt4XX/+Wj3OHTkDhIsVrc0aha+eWf0e2la33hTWZ6ouGGZyWtCoAZWVloXv37pg2bRomTZqEoUOHws/PDwBQV1eHrKws/Pzzz9i5cydyc3OxY8eOZtvz9fWFTCZDQUFBo+sFBQUmFzi3pK1tKhQKo2uK5HJ5h/yQdVS79oh9YcC+qMd+MGBfGNhiX1TV1mHdt2cBAIvG9EKov+X+g7k5ttgXYrm9L8ztl1Ytgt6xYwe+++47aDQaPProo1Cr1XB1dYWXlxcUCgUGDRqErVu3Yvbs2Th9+jTuvvvuZttzdXXFkCFDkJKSor+m0+mQkpKCmJiY1pTWoW0SEREZs/GH88grq0bXzm54YnQvscuhVmj1GqDIyEhs2bIF7777Lk6cOIFLly7hxo0b8PX1xcCBA+Hr69uq9uLj4zFnzhwMHToUUVFR2LBhAyorKzFv3jwAwOzZsxEcHIy1a9cCqF/knJWVpf/71atXcezYMXh6eqJ3795mtUlERNRe2cWV2PJjNgBg9QPhUMplIldErdHqAKTT6fDKK6/gP//5D2prazF27FgkJCSYtS7ImOnTp6OoqAirV69Gfn4+Bg4ciKSkJP0i5pycHEilhoGq3NxcDBo0SP/1+vXrsX79eowePRoHDhwwq00iIqL2evfgBdRqdbi7jx/GhfP3i71pdQD6+9//jjVr1iA2NhZubm544403UFhY2K49dpYsWYIlS5YYfa0h1DQICQkxuReRuW0SERG115oH+yOokxvujwhsskkw2b5Wb4S4Y8cOvPPOO/j222+xd+9e7Nu3D7t27eKp8ERE5FSUchn+OjYMvfw8xS6F2qDVASgnJwf33Xef/uvY2FhIJBLk5uZatDAiIiJbdOJKKeq0/I9+e9fqKbC6ujoolY3PN5HL5dyVkoiIHJJWJyA9uwSFFdWQSICnPz2Onn5e2Dk/ikde2LFWByBBEDB37txGe+ZUV1fjiSeegIeHh/5aW0+GJyIishVJmXlI3JeFvLLqRtc1Wh18PFxFqooswSJHYcyaNcsixRAREdmKpMw8LNqZAWOP3ZwvvI5vT+ZjwoBAq9dFltGmw1CJiIgcmVYnIHFfltHwAwASAIn7sjAuXA2ZlE+A2aNWL4ImIiJydOnZJU2mvW4lAMgrq0Z6don1iiKLYgAiIiK6TWGF6fDTlvvI9jAAERER3cbfS9nyTa24j2wPAxAREdFtokJ9EKhSwtTqHgmAQJUSUaE+1iyLLIgBiIiI6DYyqQTL4/oCQJMQ1PB1wqRwLoC2YwxAREREtymsqMbar0/h3n7+CPBuvNmhWqXEplmD+Qi8nWv1Y/BERESOTBAEPP3pCRRfr8XV0hv4fvkYHL9chsKKavh71U97ceTH/jEAERER3WL7Lxdx8GwRFC5SvDlzENxdXRDTq4vYZZGFcQqMiIjopjP5FfjHN6cBAM/ddwf6BHiJXBF1FAYgIiIiANUaLZbu+RW1dTrc09cPs2N6iF0SdSAGICIiIgCvfHsGp/Mr0MXDFS//KRISCdf5ODIGICIiIgARXVXwUrjg5T9FwM9L0fIbyK5xETQRERGAyQODMaaPP1TucrFLISvgCBARETktQRBQdkOj/5rhx3kwABERkdPac/gyYl87iB/PFoldClkZAxARETml34uu44V9WSiqqMHp/HKxyyErYwAiIiKnU1unw9I9x3BDo8XwXl3w55E9xS6JrIwBiIiInM6G787it6tlULnJ8eojkZDyaAunwwBEREROJe33a9h08AIAYN3DdyJQ5SZyRSQGBiAiInIaZTc0+H8fH4MgAI8M7YqJd/JEd2fFfYCIiMihaXUC0rNLUFhRDS+FC6JCffDr5VIkTOovdmkkIgYgIiJyWEmZeUjcl4W8smr9tUCVEsvH94GHgr8CnRn/6RMRkUNKyszDop0ZEG67nl9WjeWfnoCHwgUTBnAKzFlxDRARETkcrU5A4r6sJuEHgP5a4r4saHXG7iBnwABEREQOJz27pNG01+0EAHll1UjPLrFeUWRTGICIiMjhFFaYDj9tuY8cDwMQERE5HH8vpUXvI8fDAERERA7HQyFDc3s7S1D/NFhUqI+1SiIbwwBEREQORasTsOzjY/rFzrcHoYavEyaFQ8YjMJwWAxARETkUmVSCN6YPwsjevnj1kUioVY2nudQqJTbNGsxH4J0c9wEiIiKHc2dXFXb+ORoAMGVgsH4naH+v+mkvjvwQAxARETmEt1LOYWSYLwZ179zoukwqQUyvLiJVRbaKU2BERGT3dqVdwqvJZzFzyyEUlPPRdmoZAxAREdm1H88WYfUXJwEAS+7pjQBvPtpOLWMAIiIiu3UmvwKLd2VAqxPw8OBgLL6nt9glkZ1gACIiIrtUVFGD//nwMCpq6hAV6oO1D98JiYSLm8k8DEBERGR3qjVaLNhxBFdLbyDU1wPvzhoChYtM7LLIjjAAERGR3dEJAvy9FOjkLsfWucPQ2cNV7JLIzvAxeCIisjvuri7YNGsIckqqEOrrIXY5ZIcYgIiIyGZpdQLSsktwtFiCLtklULkrMSDYGxKJBDKphOGH2owBiIiIbFJSZh4S92Uhr6wagAw7zh0BAESH+mDnn6Mhl3EVB7Udf3qIiMjmJGXmYdHOjJvhp7G07BIknywQoSpyJAxARERkU7Q6AYn7svSnud9OAuDFr7Kg1Zm6g6hlDEBERGRT0rNLjI78NBAA5JVVIz27xHpFkcNhACIiIptSWGHeWV7m3kdkDAMQERHZFH8v887yMvc+ImMYgIiIyKZEhfogUKWEqUMtJAACVUpEhfpYsyxyMAxARERkMworqnG5pAoJk8IBoEkIavg6YVI4ZFKe+0VtxwBEREQ2Ib+sGjPePYSZWw6hf5AKm2YNhlrVeJpLrVJi06zBmDAgUKQqyVFwI0QiIhLd1dIbeHTLIVy6VoXgTm7QCQImDAjEuHA1Us8XYv9PaRg/Khoxvf058kMWwQBERESiyrlWhZlbDuFq6Q1093HH7gXR6NrZHQAgk0oQHeqDa6cERIf6MPyQxTAAERGRaH4vuo5Ht6Qhv7waPX09sGtBNAJVbmKXRU6AAYiIiERxoeg6Zrx3CEUVNQjz98SuBdF8tJ2shgGIiIhE0cXDFb6eCvh6KrBzfhS6eCrELomcCAMQERGJopO7K3bOj4JMKkEnd1exyyEnwwBEREQdQqsTkJ5dgsKKavh71W9cePxKKU7lleOx6B4AwFEfEg0DEBERWVxSZh4S92U1OtS0i4crrtfUoaZOB19PBeL6q0WskJwdAxAREVlUUmYeFu3MgHDb9WuVtQCAvgGeGBXma/3CiG5hEztBb9y4ESEhIVAqlYiOjkZ6enqz93/66afo168flEol7rzzTnz99deNXp87dy4kEkmjPxMmTOjIj0BERKif9krcl9Uk/NyqrLoOCheZ1WoiMkb0APTxxx8jPj4eCQkJyMjIQGRkJOLi4lBYWGj0/l9++QUzZ87E/Pnz8euvv2LKlCmYMmUKMjMzG903YcIE5OXl6f989NFH1vg4REROLT27pNG0lzH5ZdVIzy6xUkVExok+Bfbaa69hwYIFmDdvHgBg8+bN+Oqrr7B161Y8++yzTe5/4403MGHCBDz99NMAgBdffBHJycl4++23sXnzZv19CoUCarX588s1NTWoqanRf11eXg4A0Gg00Gg0bfpsxjS0Zck27RX7woB9UY/9YGCvfZFXWmn2fRqNt1n32mtfdAT2hYGpvjC3bySCIDQ3Utmhamtr4e7ujn/961+YMmWK/vqcOXNQWlqKL774osl7unfvjvj4eCxbtkx/LSEhAXv37sXx48cB1E+B7d27F66urujcuTPuvfdevPTSS+jSpYvJWtasWYPExMQm13fv3g13d/e2f0giIidyrkyCt7Nant5aEq5FmEq0Xz/kwKqqqvDoo4+irKwM3t6mQ7aoI0DFxcXQarUICAhodD0gIACnT582+p78/Hyj9+fn5+u/njBhAh5++GGEhobiwoULeO655zBx4kSkpqZCJjP+L+bKlSsRHx+v/7q8vBzdunXD+PHjm+3A1tJoNEhOTsa4ceMgl8st1q49Yl8YsC/qsR8M7LUvzhVcx8eXj6C4otboOiAJALVKgSXT7zb7XC977YuOwL4wMNUXDTM4LRF9CqwjzJgxQ//3O++8ExEREejVqxcOHDiAsWPHGn2PQqGAQtF0Pwq5XN4hP2Qd1a49Yl8YsC/qsR8M7KkvUk4VYNmeY/D1ckVxRS0kQKMQ1BB3Eib1h1LR+o0P7akvOhr7wuD2vjC3X0RdBO3r6wuZTIaCgoJG1wsKCkyu31Gr1a26HwB69uwJX19fnD9/vv1FExFRIzqdgLdSzuHPO46goqYOfp5KrJ8WAbWq8bleapUSm2YNxoQBgSJVSmQg6giQq6srhgwZgpSUFP0aIJ1Oh5SUFCxZssToe2JiYpCSktJoDVBycjJiYmJMfp8rV67g2rVrCAzkv3RERJZ0vaYOyz85jqST9csQZsf0wKr7w+HqIsWUQV2b7ARt7rQXUUcTfQosPj4ec+bMwdChQxEVFYUNGzagsrJS/1TY7NmzERwcjLVr1wIAli5ditGjR+PVV1/F/fffjz179uDIkSN47733AADXr19HYmIipk6dCrVajQsXLmDFihXo3bs34uLiRPucRESO5mJxJRb+3xGcLbgOV5kUL0zujxlR3fWvy6QSxPQy/fAJkZhED0DTp09HUVERVq9ejfz8fAwcOBBJSUn6hc45OTmQSg0zdcOHD8fu3buxatUqPPfccwgLC8PevXsxYMAAAIBMJsOJEyewfft2lJaWIigoCOPHj8eLL75odI0PERG1niAIWP7pcZwtuA5/LwU2zRqCIT06i10WkdlED0AAsGTJEpNTXgcOHGhybdq0aZg2bZrR+93c3PDtt99asjwiIrqNRCLBK9MikbjvJP45NQIB3sqW30RkQ2wiABERke25/TT3AcHeOHyxBPf2qx+hD/X1wIfzokSukqhtGICIiKgJY6e5u0gl0OoE7JgfhVFhfiJWR9R+op8FRkREtqXhNPfbz/Sq0wkQABy/XCpKXUSWxABERER65pzmvistB1odj7Eg+8YAREREeuac5p7H09zJATAAERGRXmFF8+GntfcR2SoGICIiAgDUaXXw9zLvcXZz7yOyVQxAREROTqcTsPPQJcS+dhBh/p4IVClh6sAKCYBAVf2xFkT2jAGIiMiJXS6pwqwP0rBqbyYuXqvC7vQcJEwKB4AmIchwmns4z/Qiu8d9gIiInJBOJ2Bn2iWs++Y0qmq1UMqlWBHXD3OGh0AmlWDTrMFN9gFSq5RImBTO09zJITAAERE5mUvXKrHiXyeQdvNJrqhQH7w8NQIhvh76eyYMCMS4cDVPcyeHxQBEROSAbj/G4tbwsvng70jLLoGbXIZnJ/bD43f1gNRIsOFp7uTIGICIiByMsWMs1Col1tycvnp2Qj9cr6nD0+P7onsXdxErJRIPF0ETETkQU8dY5JdVY9HODCRl5kHlLsdbMwcx/JBTYwAiInIQLR1jIQBI3JfFYyyIwABEROQweIwFkfkYgIiIHERe6Q2z7uMxFkQMQEREDiNAxWMsiMzFAEREZKcEQcCPZ4tQrdECAO7q2QU+Hq4m7+cxFkQGDEBERHbotytleOz9NMzemo7tv1wEUL9vzz8eGgAJeIwFUUsYgIiIbIxWJyAtuwRHiyVIyy5p9NTW5ZIqLN3zKya9/TN+uXANrjIpqjU6/esTBgRi06zBUN82HaZWKbFp1mAeY0F0EzdCJCKyIY03MZRhx7kjCFQp8dS4PjidX4EdqZdQq60PPA8NCkb8uD7o5tN4Px8eY0HUMgYgIiIb0bCJ4e279OSXVWP5v07ovx7RuwtWTrwDA4JVJtviMRZEzWMAIiKyAc1tYthwTS6T4L3Hh2JMXz9IJBzNIWoPrgEiIrIB5mxiqNEKUMplDD9EFsAARERkAy4UXTfrPm5iSGQZnAIjIhLR+cIKbPkxG//OuGLW/dzEkMgyGICIiCxMqxPMfgLr25MF+PjIZQD1a3w0WuMHlUpQ/yg7NzEksgwGICIiC2r8GHu9QJUSCZPCcW+/AHx5IhcB3kqM6O0LAHgsujtO5ZVj3ogQFFXUYNHODABotBiamxgSWR4DEBGRhTT3GPsTOzOgcnNB2Y06DOnRWR+AOrm74u1HB+vv3TRrcJMApb4ZoLiJIZHlMAAREVmAOY+xl92og6+nK+7t5w+tTjA6mtOwiWHq+ULs/ykN40dFI6a3P0d+iCyMAYiIyALMeYwdAF57JBJ39/Fv9h6ZVILoUB9cOyUgmjs4E3UIPgZPRGQBJ66UmnXfH1Waji2EiMzCESAioja6dr0G/zmei88yruK3q2VmvYePsRPZBgYgIqJbtPQIe02dFt+fKsS/M67iwJlC1N08qV0mARRyGW7Uao2uA+Jj7ES2hQGIiOim5h5hnzAgEFqdgHteOYDcW16P6KrC1MFdMSkyCOnZ17BoZwYk4GPsRLaOAYiICKYfYc8rq8ainRnYNGswJgwIxF09u+CXC9fw0OBgPDwoGGEBXvp7JwwI5GPsRHaCAYiInF5zj7AD9aM5ifuyMC5cjYRJ/eGpdDE5ktPwGLu5O0ETkTgYgIjI6SVl5rX4CHteWTXSs0sQ06tLi+3JpBKz7iMi8TAAEZHDMPcMLo1WhxsaLbyVcgDAydxys9rnSexEjoMBiIgcQksLmMtuaHDgTCFSThXiwJlCTBvaDc8/EA4AiOnZBe8cuNDi9+Aj7ESOgwGIiOxecwuYn9iZgT4Bnvi9qFL/yDoAHL30h/7vw3v7IlClRH5ZNR9hJ3ISDEBEZNdaWsAMAGcLrgMAwvw9ERsegNg7/DGwW2f96zKpBAmTwvkIO5ETYQAiIpug1QlIyy7B0WIJumSXmHUAaLVGi52pl8w6g2vD9EhMGdTV5Ot8hJ3IuTAAEZHoGq/fkWHHuSON1u/c6rcrZdiflY9Dv1/D8ctlqNXqzPoeEknLozd8hJ3IeTAAEZGoTK3fyb+5AeGTY3tj/sieULnVP7GVcroAb31/Xn9fZ3e5WQeMmruAmY+wEzkHBiAisghzH0G//T2m1u80XHsz5Tz6BXjhvoggAMDoPn64dK0Kd/X0QXRoF3Tt7IZRL//ABcxE1CoMQETUbi09gm6MINQHJnPW75zMK9cHoEHdO2NQ986NXucCZiJqLQYgImqXlqawNs0ajLv7+CErtxwnrpTht6tlOH6lFI8M7YZAlXnTUn1uOW/LGC5gJqLWYgAiojZNXzW8r6UprCW7f4VWJzS55/jlUkR2DTGrPnPW73ABMxG1BgMQkZNry/QVUB9+/nPsaotTWA2bDwZ4K3BncCdEdFXhzq4q3BmsQmd3V4tuQMgFzERkLgYgIjvX1tEbwLzpq7j+av0j5IIg4KlPjuNUfgUuFF1HbZ15j6C/MLk/ZseEGH2N63eISAwMQER2rK2jN4D501fDQjrjo4UxAOr30vn1cimyiysBAHKpBBpdc3sw1wvzN72Gh+t3iEgMDEBEImvLDsiAeaM3t4eHao0WOSVV+KOyFjoBZk1fHbtcCkEQ9KNAK+L6wkUmRZ8ATwSq3DD6lfY/gt6wfif1fCH2/5SG8aOize4HIqK2YAAiElFrdkC+VUujNxIAz332G07nV+ByyQ3klFTi0rUqFFbUAKjfPHDNg/3NqnF5XN9GX0+8s3FdlprCkkkliA71wbVTAqK5eJmIOhgDEFE7dPT6m9tDUG2dDgXl1fjuVEGzozcCgJIqDTZ8d67Ja15KFwR3dkMnd7lZdYYHqpo9RoJTWERkjxiAyGm1J7wAHb/+ZsW/TuDopT/wt/vD9a8t2HEEB88WmV1jdKgPRoX5onsXD/TwcUePLu7o5O6qr8FST2DxEXQisjcMQGS32jv60tbw0vB+c0ZvBEFA+Y06FF2vRlFFLYqu1+CPylr0CfBqcf1NeXUdtvyUjWWxfeChcNHX6OoiRWd3OQrKa1qsc1lsH5OPhcukEos+gcVH0InInjAAUau1d+SkoY22LPxt0J4A05app1tdr67D81+cbHb9zVOfHEfif07iWqXG6GnlL/8potkaG4y9w1+/jw4ArHmwP9Y+fCd0AjDyn99bZPExp6+IyBkxADkZMad9jLdh/sLfW9/f1gBj7tTTiStlqKiuwwuT++vXvzzzrxPYe+wqalrY+0YAUFmrRWWtVn/NW+kCXy8F/DwV8PVSoLObeetv/nzLKegAoJTLAAAyieUWH3P6ioicEQOQFbV31KOhDVuf9unINsx6+unzTFy8Vonr1VqUV2tQfkODiuo6lFdrkF9WbdbU0zsHLgAAVkzoCy9lfQCRSNBi+LlV/Lg+mDqkK7p4uOqDy62fo73rbyw5esPpKyJyNgxAVtLeUY+mbdSz1rSPOcEjcV8WxoWrjQYyjVaHypq6Fkdfnv7XCfx6uRQ6ndBo8e+a/5zEf88Xo6SyBtcqNSbrFACUVNZi3TdnTN5jjlFhvk1OHF8W2weL7+mNcwUV+J/tR1psY1iID4I7uRl9zVLrbzh6Q0TUNgxAViD2yIk50z6rvziJXn6e0Gjrr4QHeevv+eF0IVJ/L27xseu8smos2HEEeWXVqNZoUVVbhxu1WtzQaKHRCnCTS3FD0/wISkV1Hd49+DukEuC5++7QTz/ll1XjXOH1Zt97q2E9OiM8yBvebnJ4K+XwUrrA202OK3/cwD++PtXi+/8ypneTERH1zZPLgzq5WeTpKUuN4HD0hoio9WwiAG3cuBGvvPIK8vPzERkZibfeegtRUVEm7//000/x/PPP4+LFiwgLC8M///lP3HffffrXBUFAQkICtmzZgtLSUowYMQKbNm1CWFiYNT5OI+aEj7/tzYS3Ug6tIMBD4YLBt4w8fHHsKkpvaPDqt2daHH358VwRCspqUKvVoaZOh9qbf0pv1LY47VNYUYNxr/8IAAju5Ib/Pnuv/rUN353F8StlZn3eyyVVJoOKuedGje7ji4iunaDVCXCR1QegJff2xuyYHvi9uBKr9ma22Eb8+L5GQ4FWJ2Dbf7PbFV4s+fQUd0AmIhKH6AHo448/Rnx8PDZv3ozo6Ghs2LABcXFxOHPmDPz9/Zvc/8svv2DmzJlYu3YtHnjgAezevRtTpkxBRkYGBgwYAAB4+eWX8eabb2L79u0IDQ3F888/j7i4OGRlZUGpVFr186Vnl7QYPq5dr8Wj76cBAAZ374TP/jJC/9o/vj7V4uPODaMv358qQn5589+rOUq5FF5KOXw8XBtdHxriA5lEgozLpS22MeuuHgjx9YC7qwxuchncXGX6v5+4UorZWw+32MYTo5uOvgwIVgEAont2wcYfzrc5wFhy6smS62+4AzIRkXWJHoBee+01LFiwAPPmzQMAbN68GV999RW2bt2KZ599tsn9b7zxBiZMmICnn34aAPDiiy8iOTkZb7/9NjZv3gxBELBhwwasWrUKkydPBgDs2LEDAQEB2Lt3L2bMmGG9DwegsMK8QBLgrYCPhwIhXTwaXR/Txx+n8stxwowRmLj+Abgj0BsKuRSuMhlcXaRwdZHifEEFXvyq5WmfbXOjjI6aPP9AOLQ6wazHrmfd1cPkL/ARvf3aPXVkiQBjqfDC9TdERPZL1ABUW1uLo0ePYuXKlfprUqkUsbGxSE1NNfqe1NRUxMfHN7oWFxeHvXv3AgCys7ORn5+P2NhY/esqlQrR0dFITU01GYBqampQU2MYaSkvLwcAaDQaaDSmF922pIu7eV386p/uRPTNX/y3fr+XJt+BtOwSzNra8qLbcXf46du4VXQPFbb89DsKymuaCR4KDOrq1exn/dvEvnhyz3GTweNvE/tCp62DTmvkzRZsY2xfX7w1IxIvfX0a+beMjqlVCvxtYj+M7evb4j+zsX19MSZsFI5c+gOFFTXw91JgaI/OkEklrf7nPbS7N4D6NVMt1W5Kw/dsz8+aI2A/GLAvDNgXBuwLA1N9YW7fiBqAiouLodVqERAQ0Oh6QEAATp8+bfQ9+fn5Ru/Pz8/Xv95wzdQ9xqxduxaJiYlNru/fvx/u7u4tfxgTdALQyVWG0lrA8Gv+VgI6uQJFWYdgam2uJdq4Ty3B1nLpza9ubUOAAGBiQBW+Tfqmxc8zr48En12UorTW0IbKVcDDITpoLx3F15dabMIibQDAM+HAhXIJyjWAtxzo5V3Zqvc3kAG4BuDblgfJOlxycrLYJdgE9oMB+8KAfWHAvjC4vS+qqqrMep/oU2C2YuXKlY1GlsrLy9GtWzeMHz8e3t7ezbyzZfKQAjy55zgAY6MeErz0cCTi+gcYeafl2rgPwOCTBU1GTQJVSvxtYr8Wv/+t7azQCUZHTszV0MahC0X4PvUo7o0Zgrt6+Tn11JFGo0FycjLGjRsHudy8TRIdEfvBgH1hwL4wYF8YmOqLhhmclogagHx9fSGTyVBQUNDoekFBAdRqtdH3qNXqZu9v+N+CggIEBgY2umfgwIEma1EoFFAoFE2uy+Xydv+QPTCwK1xcZO1ac2KpNiZGBLd7zYocwMg+5gWm5toYEeaPsnMCRoT5O/2/yA0s8fPmCNgPBuwLA/aFAfvC4Pa+MLdfRA1Arq6uGDJkCFJSUjBlyhQAgE6nQ0pKCpYsWWL0PTExMUhJScGyZcv015KTkxETEwMACA0NhVqtRkpKij7wlJeXIy0tDYsWLerIj9MsSzzubIlFt9wzhoiIyAamwOLj4zFnzhwMHToUUVFR2LBhAyorK/VPhc2ePRvBwcFYu3YtAGDp0qUYPXo0Xn31Vdx///3Ys2cPjhw5gvfeew8AIJFIsGzZMrz00ksICwvTPwYfFBSkD1liscTjzgwwRERE7Sd6AJo+fTqKioqwevVq5OfnY+DAgUhKStIvYs7JyYFUKtXfP3z4cOzevRurVq3Cc889h7CwMOzdu1e/BxAArFixApWVlVi4cCFKS0sxcuRIJCUlWX0PICIiIrJNogcgAFiyZInJKa8DBw40uTZt2jRMmzbNZHsSiQQvvPACXnjhBUuVSERERA5E2vItRERERI6FAYiIiIicDgMQEREROR0GICIiInI6DEBERETkdBiAiIiIyOkwABEREZHTsYl9gGyRINQfOWruoWrm0mg0qKqqQnl5udOf48K+MGBf1GM/GLAvDNgXBuwLA1N90fB7u+H3uCkMQCZUVFQAALp16yZyJURERNRaFRUVUKlUJl+XCC1FJCel0+mQm5sLLy8vSCStP7PLlPLycnTr1g2XL1+Gt7e3xdq1R+wLA/ZFPfaDAfvCgH1hwL4wMNUXgiCgoqICQUFBjY7Suh1HgEyQSqXo2rVrh7Xv7e3t9D+8DdgXBuyLeuwHA/aFAfvCgH1hYKwvmhv5acBF0EREROR0GICIiIjI6TAAWZlCoUBCQgIUCoXYpYiOfWHAvqjHfjBgXxiwLwzYFwbt7QsugiYiIiKnwxEgIiIicjoMQEREROR0GICIiIjI6TAAERERkdNhABLJxYsXMX/+fISGhsLNzQ29evVCQkICamtrxS5NFH//+98xfPhwuLu7o1OnTmKXY1UbN25ESEgIlEoloqOjkZ6eLnZJovjxxx8xadIkBAUFQSKRYO/evWKXJIq1a9di2LBh8PLygr+/P6ZMmYIzZ86IXZYoNm3ahIiICP1GdzExMfjmm2/ELkt069atg0QiwbJly8QuxerWrFkDiUTS6E+/fv3a1BYDkEhOnz4NnU6Hd999FydPnsTrr7+OzZs347nnnhO7NFHU1tZi2rRpWLRokdilWNXHH3+M+Ph4JCQkICMjA5GRkYiLi0NhYaHYpVldZWUlIiMjsXHjRrFLEdXBgwexePFiHDp0CMnJydBoNBg/fjwqKyvFLs3qunbtinXr1uHo0aM4cuQI7r33XkyePBknT54UuzTRHD58GO+++y4iIiLELkU0/fv3R15env7Pzz//3LaGBLIZL7/8shAaGip2GaLatm2boFKpxC7DaqKiooTFixfrv9ZqtUJQUJCwdu1aEasSHwDh888/F7sMm1BYWCgAEA4ePCh2KTahc+fOwvvvvy92GaKoqKgQwsLChOTkZGH06NHC0qVLxS7J6hISEoTIyEiLtMURIBtSVlYGHx8fscsgK6mtrcXRo0cRGxurvyaVShEbG4vU1FQRKyNbUlZWBgBO//8NWq0We/bsQWVlJWJiYsQuRxSLFy/G/fff3+j/M5zRuXPnEBQUhJ49e+Kxxx5DTk5Om9rhYag24vz583jrrbewfv16sUshKykuLoZWq0VAQECj6wEBATh9+rRIVZEt0el0WLZsGUaMGIEBAwaIXY4ofvvtN8TExKC6uhqenp74/PPPER4eLnZZVrdnzx5kZGTg8OHDYpciqujoaHz44Yfo27cv8vLykJiYiFGjRiEzMxNeXl6taosjQBb27LPPNlmgdfuf23+5Xb16FRMmTMC0adOwYMECkSq3vLb0BREZLF68GJmZmdizZ4/YpYimb9++OHbsGNLS0rBo0SLMmTMHWVlZYpdlVZcvX8bSpUuxa9cuKJVKscsR1cSJEzFt2jREREQgLi4OX3/9NUpLS/HJJ5+0ui2OAFnYU089hblz5zZ7T8+ePfV/z83NxT333IPhw4fjvffe6+DqrKu1feFsfH19IZPJUFBQ0Oh6QUEB1Gq1SFWRrViyZAm+/PJL/Pjjj+jatavY5YjG1dUVvXv3BgAMGTIEhw8fxhtvvIF3331X5Mqs5+jRoygsLMTgwYP117RaLX788Ue8/fbbqKmpgUwmE7FC8XTq1Al9+vTB+fPnW/1eBiAL8/Pzg5+fn1n3Xr16Fffccw+GDBmCbdu2QSp1rAG51vSFM3J1dcWQIUOQkpKCKVOmAKif8khJScGSJUvELY5EIwgCnnzySXz++ec4cOAAQkNDxS7Jpuh0OtTU1IhdhlWNHTsWv/32W6Nr8+bNQ79+/fDMM884bfgBgOvXr+PChQt4/PHHW/1eBiCRXL16FWPGjEGPHj2wfv16FBUV6V9zxv/6z8nJQUlJCXJycqDVanHs2DEAQO/eveHp6SlucR0oPj4ec+bMwdChQxEVFYUNGzagsrIS8+bNE7s0q7t+/Xqj/4rLzs7GsWPH4OPjg+7du4tYmXUtXrwYu3fvxhdffAEvLy/k5+cDAFQqFdzc3ESuzrpWrlyJiRMnonv37qioqMDu3btx4MABfPvtt2KXZlVeXl5N1oB5eHigS5cuTrc2bPny5Zg0aRJ69OiB3NxcJCQkQCaTYebMma1vzCLPklGrbdu2TQBg9I8zmjNnjtG++OGHH8QurcO99dZbQvfu3QVXV1chKipKOHTokNglieKHH34w+jMwZ84csUuzKlP/v7Bt2zaxS7O6//mf/xF69OghuLq6Cn5+fsLYsWOF/fv3i12WTXDWx+CnT58uBAYGCq6urkJwcLAwffp04fz5821qSyIIgtCeNEZERERkbxxr0QkRERGRGRiAiIiIyOkwABEREZHTYQAiIiIip8MARERERE6HAYiIiIicDgMQEREROR0GICIiInI6DEBE1KHGjBmDZcuWteo9c+fOhUQigUQiwd69ezukLlvhTJ+VyJYwABGRTZowYQLy8vIwceJEq37fkJAQSCQSHDp0qNH1ZcuWYcyYMU3uv3LlClxdXU2eydQQbiQSCVQqFUaMGIHvv/9e//obb7yBvLw8i34GImoZAxAR2SSFQgG1Wg2FQmHRdufOnYs1a9Y0e49SqcQzzzxjVnsffvghHnnkEZSXlyMtLc3oPdu2bUNeXh7++9//wtfXFw888AB+//13APWHnDrjAchEYmMAIqIOp9PpsGLFCvj4+ECtVrcYQIy5ePEiJBIJPvnkE4waNQpubm4YNmwYzp49i8OHD2Po0KHw9PTExIkTUVRU1K56Fy5ciEOHDuHrr79u9j5BELBt2zY8/vjjePTRR/HBBx8Yva9Tp05Qq9UYMGAANm3ahBs3biA5ObldNRJR+zAAEVGH2759Ozw8PJCWloaXX34ZL7zwQpsDQEJCAlatWoWMjAy4uLjg0UcfxYoVK/DGG2/gp59+wvnz57F69ep21RsaGoonnngCK1euhE6nM3nfDz/8gKqqKsTGxmLWrFnYs2cPKisrm23bzc0NAFBbW9uuGomofRiAiKjDRUREICEhAWFhYZg9ezaGDh2KlJSUNrW1fPlyxMXF4Y477sDSpUtx9OhRPP/88xgxYgQGDRqE+fPn44cffmh3zatWrUJ2djZ27dpl8p4PPvgAM2bMgEwmw4ABA9CzZ098+umnJu+vqqrCqlWrIJPJMHr06HbXSERtxwBERB0uIiKi0deBgYEoLCxsd1sBAQEAgDvvvLPRtVvb3rVrFzw9PfV/du3ahX/84x+Nrv30009Nvo+fnx+WL1+O1atXGx2tKS0txWeffYZZs2bpr82aNcvoNNjMmTPh6ekJLy8v/Pvf/8YHH3zQpE+IyLpcxC6AiByfXC5v9LVEIml2asnctiQSidFrt7b94IMPIjo6Wv/1M888g+DgYPz1r3/VXwsODjb6veLj4/HOO+/gnXfeafLa7t27UV1d3ahtQRCg0+lw9uxZ9OnTR3/99ddfR2xsLFQqFfz8/FrzcYmogzAAEZFD8/LygpeXV6OvfXx80Lt37xbf6+npieeffx5r1qzBgw8+2Oi1Dz74AE899RTmzp3b6Ppf/vIXbN26FevWrdNfU6vVZn0/IrIeToERETVj4cKFUKlU2L17t/7asWPHkJGRgT//+c8YMGBAoz8zZ87E9u3bUVdXJ2LVRNQSBiAiombI5XK8+OKLqK6u1l/74IMPEB4ejn79+jW5/6GHHkJhYWGLj9ATkbgkgiAIYhdBRHSruXPnorS01KmOhpBIJPj8888xZcoUsUshcgocASIim/Tll1/C09MTX375pdildKgnnngCnp6eYpdB5HQ4AkRENqewsBDl5eUA6h+Z9/DwELmijuNMn5XIljAAERERkdPhFBgRERE5HQYgIiIicjoMQEREROR0GICIiIjI6TAAERERkdNhACIiIiKnwwBERERETocBiIiIiJzO/weSZTC8pntdCgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(h, fc_pf, 'o--')\n",
    "plt.grid()\n",
    "plt.xlabel(\"h [m+NAP]\")\n",
    "plt.ylabel(\"P(Z<0|h)\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0dfcd735",
   "metadata": {},
   "source": [
    "The load $h$ is the integrand and follows an exponential distribution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "id": "eccba97a",
   "metadata": {},
   "outputs": [],
   "source": [
    "integrand = Stochast()\n",
    "integrand.name = \"h\"\n",
    "integrand.distribution = DistributionType.exponential\n",
    "integrand.shift = 0.5\n",
    "integrand.scale = 1.0"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "74c10f0c",
   "metadata": {},
   "source": [
    "The integration of the fragility curve with the water level statistics is performed as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "id": "2a5573cb",
   "metadata": {},
   "outputs": [],
   "source": [
    "dp = fragility_curve.integrate(integrand)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d50177bd",
   "metadata": {},
   "source": [
    "The results are:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "id": "ac71a738",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Reliability:\n",
      " Reliability index = 1.8853\n",
      " Probability of failure = 0.0297\n",
      " Convergence = 0.0 (converged)\n",
      " Model runs = 16001\n",
      "Alpha values:\n",
      " h: alpha = -0.5209, x = 2.3137\n",
      " t_p: alpha = -0.8555, x = 9.6084\n",
      " h_s: alpha = -0.4226, x = 3.6859\n",
      " h_crest: alpha = 0.0168, x = 9.9983\n",
      "\n"
     ]
    }
   ],
   "source": [
    "dp.print()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b836b711",
   "metadata": {},
   "source": [
    "Reliability calculations performed without the fragility curve, but instead using the distribution function for the water level $h$, should yield results that are close to those obtained using the fragility curve.\n",
    "\n",
    "We begin with the `form` method, which provides a reliability index that deviates from the one obtained from the integration of the water level distribution with the fragility curve. It is expected that the method has converged to a local optimum in this case."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "id": "f98511d9",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Reliability (FORM)\n",
      " Reliability index = 2.0502\n",
      " Probability of failure = 0.0202\n",
      " Convergence = 0.0089 (converged)\n",
      " Model runs = 165\n",
      "Alpha values:\n",
      " t_p: alpha = -0.7813, x = 9.5736\n",
      " tan_alpha: alpha = 0.0, x = 0.3333\n",
      " h_s: alpha = -0.3859, x = 3.6793\n",
      " h_crest: alpha = 0.0154, x = 9.9983\n",
      " h: alpha = -0.4904, x = 2.3491\n",
      "\n"
     ]
    }
   ],
   "source": [
    "project.variables[\"h\"].distribution = DistributionType.exponential\n",
    "project.variables[\"h\"].shift = 0.5\n",
    "project.variables[\"h\"].scale = 1.0\n",
    "\n",
    "project.run()\n",
    "\n",
    "project.design_point.print()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4e441dce",
   "metadata": {},
   "source": [
    "The results are improved, when the `crude_monte_carlo` method is used:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "id": "2e0976b8",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Reliability:\n",
      " Reliability index = 1.8832\n",
      " Probability of failure = 0.0298\n",
      " Convergence = 0.0255 (not converged)\n",
      " Model runs = 50001\n",
      "Alpha values:\n",
      " t_p: alpha = -0.7713, x = 9.1209\n",
      " tan_alpha: alpha = 0.0, x = 0.3333\n",
      " h_s: alpha = -0.3697, x = 3.5677\n",
      " h_crest: alpha = 0.0401, x = 9.9961\n",
      " h: alpha = -0.5164, x = 2.2994\n",
      "\n"
     ]
    }
   ],
   "source": [
    "project.settings.reliability_method = ReliabilityMethod.crude_monte_carlo\n",
    "project.settings.minimum_samples = 1000\n",
    "project.settings.maximum_samples = 50000\n",
    "project.settings.variation_coefficient = 0.02\n",
    "\n",
    "project.run()\n",
    "\n",
    "project.design_point.print()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "base",
   "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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}