six_unit_weight_sampling.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 117,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from scipy import stats\n",
    "import numpy as np\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [],
   "source": [
    "f_mean = 22"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "22\n"
     ]
    }
   ],
   "source": [
    "print(f_mean)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [],
   "source": [
    "f_cov = 0.10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [],
   "source": [
    "f_std = f_mean * f_cov"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.2\n"
     ]
    }
   ],
   "source": [
    "print(f_std)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Monte Carlo simulations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {},
   "outputs": [],
   "source": [
    "simulations_number = 100_000 #assumed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.4359949 , 0.02592623, 0.54966248, 0.43532239, 0.4203678 ])"
      ]
     },
     "execution_count": 79,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "u = np.random.uniform(low=0, high=1, size=simulations_number) #probability realisations\n",
    "u[0:5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-0.16113153, -1.9443566 ,  0.12480864, -0.16283953, -0.20095264])"
      ]
     },
     "execution_count": 80,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = stats.norm.ppf(u) #inverse of the normal distribution\n",
    "x[0:5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([21.64551063, 17.72241548, 22.27457902, 21.64175304, 21.55790419])"
      ]
     },
     "execution_count": 81,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "r = x*f_std + f_mean #vector of realisations\n",
    "r[0:5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [],
   "source": [
    "compressive_resistance_distribution = stats.norm.pdf(np.sort(r), f_mean, f_std)\n",
    "#https://stackoverflow.com/questions/10138085/how-to-plot-normal-distribution\n",
    "#https://www.askpython.com/python/normal-distribution\n",
    "#https://vitalflux.com/normal-distribution-explained-python-examples/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {},
   "outputs": [],
   "source": [
    "# x_label = input()\n",
    "# y_label = input()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 103,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Probability density')"
      ]
     },
     "execution_count": 103,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEMCAYAAADTfFGvAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAA5/UlEQVR4nO3deXxU9bn48c+TyQZhCSEhgbAEMCyRnYBY3BUErOJSrbZVa73X0mqttt5ee7v5u729tYt6bWvrUrXaRat1Q8UK7ooo+xbWAAECIQlrEiDrPL8/zokdxywzkJkzyTzv1+u85sw553vmmcnJPHO+53u+X1FVjDHGmFAleB2AMcaYzsUShzHGmLBY4jDGGBMWSxzGGGPCYonDGGNMWCxxGGOMCUui1wFEQ2Zmpubl5XkdhjHGdCorVqzYr6pZwcvjInHk5eWxfPlyr8MwxphORUR2trTcqqqMMcaExRKHMcaYsFjiMMYYExZLHMYYY8JiicMYY0xYLHEYY4wJS1w0xzUmVqgq2/cfZc+h4/gShCF9uzOwT3evwzImLJY4jImCI8caeGJJCX/7eBf7qmo/tW5YZhrXnT6Eq6cOJjXJ51GExoTOEocxEbawaB//9cJ69tfUcc7ILL59QT7Ds3rgV2VTWRUvry3jrpc38OSSnfzqynFMHpLhdcjGtEniYQTAwsJCtTvHTbSpKr95s5j73thCQf9e/OKKcYwd2LvFbd/fWsl/vbCOssO1/O9lY7lqyqAoR2vMZ4nIClUtDF5uF8eNiQBV5b9f2cB9b2zh8km5vHjz9FaTBsCZ+Vm88q0zOX14X7733Fr+tHhHFKM1JjwRTRwiMktENotIsYjc2cL6USKyRETqROSOgOUjRWR1wFQlIre56+4SkT0B6+ZE8j0YcyJ++1Yxjy8u4WvTh3LPleNJTmz/X613tyQe++oUZhRkc9fLG3huRWkUIjUmfBFLHCLiAx4AZgMFwDUiUhC02UHgVuDXgQtVdbOqTlDVCcBk4BjwQsAm9zWvV9UFkXoPxpyIBevKuHfRFq6YNJAfXjQaEQm5bJIvgd99aSKfG96X7z+/jtW7D0cuUGNOUCTPOKYCxaq6XVXrgaeBuYEbqGqFqi4DGtrYz/nANlVtsZdGY2LJ9soavvePtUwcnM7PLx9LQkLoSaNZSqKPB740iX69Upj35xUcqKmLQKTGnLhIJo5cYHfA81J3WbiuBp4KWnaLiKwVkcdEpE9LhUTkJhFZLiLLKysrT+BljQlPQ5Ofbz21iiSf8MCXJoVUPdWaPmnJPPiVyRw4WsePXyrqwCiNOXmRTBwt/dQKqwmXiCQDlwDPBiz+AzAcmACUAfe0VFZVH1bVQlUtzMr6zDgkxnS4h9/bTtHeKn5++VgGpHc76f2Nye3N7TNG8Oq6Ml5es7cDIjSmY0QycZQCgW0KBwLhHv2zgZWqWt68QFXLVbVJVf3AIzhVYsZ4qriihvvf3MqcsTnMGtO/w/Z705nDmDg4nR+/tJ7Dx+o7bL/GnIxIJo5lQL6IDHXPHK4G5oe5j2sIqqYSkcD/ysuA9ScVpTEnSVX5yfz1dEvycdclp3bovhN9CfzvZWM5cryBexdt6dB9G3OiIpY4VLURuAV4HdgIPKOqRSIyT0TmAYhIjoiUAt8BfigipSLSy13XHZgBPB+061+KyDoRWQucC9weqfdgTCje3lzB4uID3H5BPv16pnb4/kf378W104bwl492srGsqsP3b0y47M5xY05CQ5OfWf/3Hqrw+u1nkeSLzG+xw8fqOffX7zC6fy/+9u/TIvIaxgSzO8eNiYCnl+5iW+VR/mvO6IglDYD07sl867x8Ptx2gA+L90fsdYwJhSUOY05QbUMTv3u7mKl5GZw/ul/EX+9Lpw2mf+9Ufr1wM/FQU2BilyUOY07Qs8t3U15Vx20X5Id1d/iJSk3y8a3z8lm56zBvb66I+OsZ0xpLHMacgLrGJn7/zjYKh/Th9OF9o/a6VxYOZHBGd+5/Y6uddRjPWOIw5gT8Y0UpZUdq+XaUzjaaJfkS+PrZw1hTeoSPth+M2usaE8gShzFhavIrD727nQmD0jnjlMyov/4VkwaS2SOZB9/dFvXXNgYscRgTtkUbytl18Bg3nTUsqmcbzVKTfNwwfSjvbqm0+zqMJyxxGBOmxz7YQW56N2YWZHsWw1dOG0Jaso+H7KzDeMAShzFhWFd6hKUlB7lheh6JEbxvoz29uyfxxSmDeXVdGRXVtZ7FYeKTJQ5jwvDoB9tJS/bFxJjg154+hIYm5emlu9vf2JgOZInDmBBVVNfyytoyriwcRK/UJK/DYWhmGmfmZ/K3j3fR2OT3OhwTRyxxGBOiZ5eX0uhXrjt9iNehfOK60/PYV1XLog3l7W9sTAexxGFMCPx+5ellu5g2LINhWT28DucT543qR256N55cYiMrm+ixxGFMCBZv28/ug8e5Zupgr0P5FF+C8OVpg1my/QDFFdVeh2PihCUOY0Lw1NJd9OmexIWn5ngdymdcOXkQvgTh2eWlXodi4oQlDmPaUVldx8Kicq6YNJDUJJ/X4XxGVs8UzhvVj+dW7qHBLpKbKLDEYUw7nlvpXBS/OsaqqQJdVTiI/TV1vLO50utQTBywxGFMG1SVZ5btZkpeH07pFzsXxYOdMzKLzB4pPLvc7ukwkRfRxCEis0Rks4gUi8idLawfJSJLRKRORO4IWlfiji2+WkSWByzPEJFFIrLVfewTyfdg4tua0iNs33+UKyYN9DqUNiX5Erh8Ui5vbaqgsrrO63BMFxexxCEiPuABYDZQAFwjIgVBmx0EbgV+3cpuzlXVCUFj3t4JvKmq+cCb7nNjIuKFlaUkJyYwe2x/r0Np15WTB9LoV15ctcfrUEwXF8kzjqlAsapuV9V64GlgbuAGqlqhqsuAhjD2Oxd4wp1/Ari0A2I15jMamvy8vLaMGaOz6d3N+zvF25Of3ZOJg9P5xwprXWUiK5KJIxcIrHAtdZeFSoGFIrJCRG4KWJ6tqmUA7mPkB3s2cendzZUcPFrPZRPDOWy9dfnEXDaXV7Npn3W3biInkomjpYEKwhnrcrqqTsKp6rpZRM4K68VFbhKR5SKyvLLSWpqY8L2wag8ZacmcPTLL61BCNmdsf3wJwour9nodiunCIpk4SoHALkQHAiEfzaq6132sAF7AqfoCKBeR/gDuY0Ur5R9W1UJVLczK6jz/+CY2HDnewKKN5Vw8rj9JHnafHq6+PVI4Mz+Tl9fsxe+3MclNZETyP2IZkC8iQ0UkGbgamB9KQRFJE5GezfPATGC9u3o+cL07fz3wUodGbQzw2roy6hv9XBbjralaMnfCAPYcPs6KXYe8DsV0UYmR2rGqNorILcDrgA94TFWLRGSeu/5BEckBlgO9AL+I3IbTAisTeMEdljMR+Juq/tPd9d3AMyJyI7ALuDJS78HErxdW7WFYZhrjB/b2OpSwzSzIITVpHS+t3sOUvAyvwzFdUMQSB4CqLgAWBC17MGB+H04VVrAqYHwr+zwAnN+BYRrzKeVVtSwtOci3z8/3ZEzxk5WWksiMghxeXVvGTy4+tVNVtZnOwY4oY4K8tq4MVbioE9y70Zq54wdw6FgD72+1hiGm41niMCbIq+vKGJndk/zsnl6HcsLOGpFF725J1rrKRIQlDmMC7DtSy7KSQ1w0rvOebQAkJyYwZ2wOb24sp7ahyetwTBdjicOYAAvWlQF0+sQBMHtMf47WN/HuFquuMh3LEocxAV5dV8bo/r0YHkPDw56o04f3Jb17Eq+5ydCYjmKJwxjX3sPHWbHzEJ/vAmcb4PSYO7Mgmzc2VlDXaNVVpuNY4jDG1VxNNacTt6YKNntsf2rqGvlg636vQzFdiCUOY1yvrivj1AG9GJqZ5nUoHWb68Ex6pSbyqlVXmQ5kicMYoPTQMVbtOtwlLooHSk5MYEZBDos2lFPfaOORm45hicMY4J/r9wGd+6a/1swZm0N1bSOLt1l1lekYljiMARZuKGdUTk+G9O061VTNzsjPpGdKIgvWWnWV6RiWOEzcO3i0nuUlB5lZkO11KBGRkujjgoJsFm4op6HJqqvMybPEYeLeW5sq8CvMKMjxOpSImT0mhyPHG1iy7YDXoZguwBKHiXsLi/bRv3cqY3J7eR1KxJw1IovuyT4WbtjndSimC2g3cYjImGgEYowXahuaeH/rfmYUZHfKLtRDlZrk4+wRWSzaUG4jA5qTFsoZx4MislREviki6ZEOyJho+mDrfo43NDGji17fCDTz1GzKq+pYu+eI16GYTq7dxKGqZwBfxhk/fLmI/E1EZkQ8MmOiYNGGcnqmJHLa0L5ehxJx543MxpcgLCyy6ipzckK6xqGqW4EfAv8JnA38RkQ2icjlkQzOmEhq8itvbCznnFH9SE7s+pf7endPYtqwDBZuKPc6FNPJhXKNY5yI3AdsBM4DLlbV0e78fRGOz5iIWbXrEAeO1nfZZrgtmVmQQ3FFDdsqa7wOxXRiofzM+h2wEhivqjer6koAVd2LcxbSKhGZJSKbRaRYRO5sYf0oEVkiInUickfA8kEi8raIbBSRIhH5dsC6u0Rkj4isdqc5ob5ZYwIt2lBOkk84Z2SW16FETfO1nEV21mFOQiiJ43lV/bOqHm9e0PxFrqp/bq2QiPiAB4DZQAFwjYgUBG12ELgV+HXQ8kbgu+6ZzTTg5qCy96nqBHdaEMJ7MOZTVJWFG8qZNqwvPVOTvA4nagakd2Nsbm9et+sc5iSEkjiua2HZV0MoNxUoVtXtqloPPA3MDdxAVStUdRnQELS8LODMphqnmiw3hNc0JiTbKmvYsf8oM0/tujf9tWZmQTardh2moqrW61BMJ9Vq4hCRa0TkZWCoiMwPmN4GQrn9NBfYHfC8lBP48heRPGAi8HHA4ltEZK2IPCYifVopd5OILBeR5ZWVNnSm+bTmC8QzRsfP9Y1mzcly0UarrjInpq0zjg+Be4BN7mPz9F1gVgj7buluqrDuPBKRHsBzwG2qWuUu/gMwHJgAlLkxffaFVB9W1UJVLczKip86bBOaRRvKGTewNzm9U70OJepGZPdgSN/uLCyyxGFOTGJrK1R1J7ATOP0E912Kc+9Hs4HA3lALi0gSTtL4q6o+HxBXecA2jwCvnGB8Jk5VVNWyatdhvjtjhNeheEJEmFmQzZ8+LKG6tiGurvGYjtFWVdUH7mO1iFQFTNUiUtVauQDLgHwRGSoiycDVwPxQghKn74dHgY2qem/QusABEy4D1oeyT2OavbGxAiAur280u/DUHBqalHc2WzWuCV9bZxxnuI89T2THqtooIrcArwM+4DFVLRKRee76B0UkB1gO9AL8InIbTgusccC1wDoRWe3u8r/cFlS/FJEJONVeJcDXTyQ+E78WbdjH4IzujMju4XUonpk4uA+ZPZJZuKGci8cP8Doc08m0mjiaichwoFRV60TkHJwv9SdV9XB7Zd0v+gVByx4MmN+HU4UV7ANavkaCql7b3usa05qaukYWbzvAtdOGdOlODdvjSxAuGJ3NK2vLqG/0x8Wd86bjhHK0PAc0icgpONVHQ4G/RTQqYyLkvS2V1Df64+pu8dbMKMimpq6Rj7bbGB0mPKEkDr+qNuJcT/g/Vb0d6HoDM5u4sGhDOX26JzF5SIutuOPK9FMy6ZZkY3SY8IWSOBpE5Brgev7VgsmaYZhOp6HJz1ubKjhvVDaJPquaaR6j440NFajaGB0mdKH899yA0yT3Z6q6Q0SGAn+JbFjGdLxlJQc5crwhLsbeCNWMgmz2VdWyzsboMGFo9+K4qm7A6U+q+fkO4O5IBmVMJCwsKiclMYGzRmR6HUrMOG9UP3eMjnLGDUz3OhzTSYTSrfp0EVkkIltEZLuI7BCR7dEIzpiOoqos2lDOmfmZdE9u9/dS3OiTlsyUvD7WW64JSyhVVY8C9wJnAFOAQvfRmE5jY1k1ew4ft2qqFswoyGFzeTU7Dxz1OhTTSYSSOI6o6mtuT7YHmqeIR2ZMB1q4YR8icN4oSxzBZtoYHSZMoSSOt0XkVyJyuohMap4iHpkxHWjRhnImD+5DVs8Ur0OJOYMyujMqp6cNKWtCFkpl72nuY2HAMsUZOtaYmLfn8HGK9lbx/dmjvA4lZs0syOZ3bxdz8Gg9GWnJXodjYly7Zxyqem4LkyUN02m80Tz2hl3faNWMghz8Cm/aGB0mBKG0qsoWkUdF5DX3eYGI3Bj50IzpGAs37OOUfj0YlhW/nRq2Z0xuL/r3TrXrHCYkoVzj+BNOD7fNXWhuAW6LUDzGdKgjxxv4ePtBO9toh4gwoyCb97ZWcry+yetwTIwLJXFkquozgB+c7tIBO7JMp/DO5goa/WqJIwQzCrKpbfDzQfF+r0MxMS6UxHFURPriDvsqItMA65/AdAoLN5ST1TOFCXZXdLtOG9qXnimJLLJOD007QmlV9R2ckfuGi8hiIAv4QkSjMqYD1DU28c6mCi6ZkEtCQvyOvRGq5MQEzh3Vjzc3VtDkV3z2mZlWhNKqaiVwNvA5nNH2TlXVtZEOzJiTtWTbAY7WN9nYG2GYUZDNgaP1rNx1yOtQTAxr9YxDRC5vZdUIEUFVn49QTMZ0iEUbyume7OP04X29DqXTOGdkFkk+YdGGcqbkZXgdjolRbZ1xXOxON+L0V/Vld/oj8JVQdi4is0Rks4gUi8idLawfJSJLRKRORO4IpayIZLidLm51H21EHvMZfr/TqeHZI7JITfJ5HU6n0TM1idOHZ7KwaJ+N0WFa1WriUNUbVPUGnIviBap6hapeAZwayo5FxAc8AMwGCoBrRKQgaLODOF22/zqMsncCb6pqPvCm+9yYT1m75wgV1XXMPNWqqcI1oyCbkgPHKK6o8ToUE6NCaVWVp6plAc/LgREhlJsKFKvqdlWtB54G5gZu4HacuAxoCKPsXOAJd/4J4NIQYjFxZtGGffgShHNH9vM6lE5nxmgn2VrfVaY1oSSOd0TkdRH5qohcD7wKvB1CuVxgd8DzUndZKNoqm92cyNzHFr8ZROQmEVkuIssrKytDfFnTVSwsKue0oRmkd7d+l8KV0zuV8QN7213kplWhtKq6BXgQGA9MAB5W1W+FsO+W2vKFWml6MmWdjVUfVtVCVS3MysoKp6jp5LZX1rC1osZaU52EGQXZrN59mPKqWq9DMTEolDMOVPUFVb3dnV4Icd+lwKCA5wOBvR1QtlxE+gO4jxUh7tPEieZfyjNOzfE4ks5rRoHz2b1hnR6aFoSUOE7QMiBfRIaKSDJwNc6NhCdbdj5wvTt/PfBSB8ZsuoBFG8oZk9uL3PRuXofSaY3I7sHgjO5WXWVaFLHE4fZpdQtOB4kbgWdUtUhE5onIPAARyRGRUpy7038oIqUi0qu1su6u7wZmiMhWYIb73BgAKqvrWLHrEDML7GzjZIgIMwuy+bD4ADV1jV6HY2JMu12OiMjngQWq6g9356q6AFgQtOzBgPl9ONVQIZV1lx8Azg83FhMf3txYjirWDLcDzCjI5o8f7ODdzZVcNK6/1+GYGBLKGcfVwFYR+aWIjI50QMacjIUbyhmU0Y2R2T29DqXTmzykD326J1mnh+YzQmlV9RVgIrANeNy90/smEbH/TBNTauoa+aB4PzMLchCxDvpOVqIvgfNHZ/PWpgoamsKucDBdWKitqqqA53BuxOsPXAasFJFQmuUaExXvbamkvtFvzXA70IyCbKpqG1m646DXoZgYEsrQsZeIyAvAW0ASMFVVZ+Pc13FHm4WNiaKFRfvISEtm8hDrvqyjnJmfSUpigrWuMp8SyhnHF4D7VHWcqv5KVSsAVPUY8LWIRmdMiBqa/Ly1qYLzR/Uj0RfJVubxpXtyImfmZ1mnh+ZTQvkPK1PV9wIXiMgvAFT1zYhEZUyYlu44SFVtIzPtpr8ON7Mgm71HainaW+V1KCZGhJI4ZrSwbHZHB2LMyVhYtI9uST7OzM/0OpQu57zR/RDBqqvMJ1pNHCLyDRFZB4wSkbUB0w7ARgA0MUNVWbihnDPzM23sjQjI7JFC4ZA+1luu+URbZxx/wxnI6SX+NajTxcBkt4muMTFh/Z4qyo7UWjVVBM0oyGZjWRW7Dx7zOhQTA9pKHKqqJcDNQHXAhIjYmJImZrxe5Iy9cf4oG3sjUqzTQxOovTMOgBXAcvdxRcBzYzynqixYX8a0YRn0SbOxNyJlaGYap/TrwcIiSxymjb6qVPXz7uPQ6IVjTHi2VtSwvfIoX5tuh2mkzSzI5qH3tnP4WL0NkBXn2ro4PqmtKZpBGtOaBevKEIEL7fpGxM0oyKbJr7y92YbAiXdt9Y57TxvrFDivg2MxJmyvrdvHlLwMsnqmeB1Klzd+YDr9eqawsKicyya22Km1iRNtVVWdG81AjAlXcUUNm8urueviAq9DiQsJCcKMgmyeX7mH4/VNdEu2ps/xqq2qqvPcx8tbmqIXojEt++f6MgBmjbGxIqLlorH9Od7QZNVVca6tqqqzcTo2vLiFdQo8H5GIjAnRa+v3MWlwOjm9U70OJW5MHZpBZo9kXl1XxpyxlrDjVVtVVT9xH2+IXjjGhGbngaMU7a3ihxfZ2GLRlOhL4MJTc3h+5R6O1TfSPbndQURNFxRKt+p9ReQ3IrJSRFaIyP0i0jeUnYvILBHZLCLFInJnC+vF3Xex253JJHf5SBFZHTBVicht7rq7RGRPwLo5Yb5n0wW8tt4ZlW7WGGtNFW0XjXOrqzZVeh2K8UgonRw+DVQCV+B0sV4J/L29QiLiAx7A6RCxALhGRIKvYs4G8t3pJuAPAKq6WVUnqOoEYDJwDHghoNx9zevdsclNnHltXRnjB/ZmYJ/uXocSd04b2pfMHsksWFfmdSjGI6EkjgxV/amq7nCn/wHSQyg3FShW1e2qWo+TgOYGbTMXeFIdHwHpIhJccXo+sE1Vd4bwmiYOlB46xprSI3ZR3CO+BGHWmBze3FTOsfpGr8MxHgglcbwtIleLSII7XQW8GkK5XGB3wPNSd1m421wNPBW07Ba3ausxEWlxuDd3XPTlIrK8stJOqbuS19Y51VSzrZrKMxeNHUBtg9+qq+JUW81xq0WkCvg6Tr9V9e70NHB7CPuWFpYFDyHW5jYikgxcAjwbsP4PwHBgAlBGKzcqqurDqlqoqoVZWVkhhGs6i5fX7mVMbi/yMtO8DiVuOa2rUnh13V6vQzEeaDVxqGpPVe3lPiaoaqI7JahqrxD2XQoMCng+EAg+ytrbZjawUlU/6VlNVctVtUlV/cAjOFViJk7s2H+UtaVHuGT8AK9DiWu+BGH2mBze2lRh1VVxKKTBmUWkj4hMFZGzmqcQii0D8kVkqHvmcDUwP2ib+cB1buuqacARVQ284nYNQdVUQddALgPWh/IeTNfw8hrnd8Xnx1ni8NpF4/pT2+CM9W7iS7uNsEXk34Bv45wNrAamAUtop68qVW0UkVuA1wEf8JiqFonIPHf9g8ACYA5QjNNy6pN7RkSkO86wtV8P2vUvRWQCTpVWSQvrTRelqsxfs5epeRkMSO/mdThxb0qeW121tswSeZwJ5e6dbwNTgI9U9VwRGQX8v1B27jaVXRC07MGAecUZKKqlsseAz9wvoqrXhvLapuvZtK+a4ooafnrpGK9DMTjVVReNzeHpZbuprm2gZ2qS1yGZKAmlqqpWVWsBRCRFVTcBIyMbljGfNX/NXnwJwhxrTRUzLpmQS12jn3+6N2Sa+BBK4igVkXTgRWCRiLzEZy9yGxNRqsrLa/Yy/ZRM+vawLtRjxaTB6QzO6M5Lq+0rIZ60mzhU9TJVPayqdwE/Ah4FLo1wXMZ8yspdhyk9dNxaU8UYEeHSCQNYvG0/5VW1XodjoiTUVlWTRORWYBxQ6t4JbkzUvLxmL8mJCVx4arbXoZggcyfmovqvFm+m6wulk8MfA0/gXKjOBB4XkR9GOjBjmjU2+XllbRnnjexnF2Bj0PCsHowb2JsXVu3xOhQTJaGccVwDTFHVn7hdrU8DvhzZsIz5l/eL97O/po5LJ1o1Vay6dEIuRXur2Fpe7XUoJgpCSRwlQOBIOSnAtohEY0wLnltRSnr3JM4d1c/rUEwrPj++PwkCL662s4540FZfVb8Vkd8AdUCRiPxJRB7HuVO7JloBmvh25HgDCzeUc8n4AaQk2hjXsapfz1TOyM/ipdV78fuDu6QzXU1bNwAudx9X8OmxMN6JWDTGBHl1bRn1jX6umDTQ61BMOy6dMIDvPLOGFbsOMSUvw+twTAS1NXTsE83zbl9TI9ynm1W1IdKBGQPw3MpSTunnXHw1se3CU3Ponrye51aUWuLo4kJpVXUOsBVnNL/fA1tC7OTQmJOyY/9RVuw8xBWTBiLSUg/8JpakpSRy0dj+vLxmL0frrMfcriyUi+P3ADNV9WxVPQu4ELgvsmEZA8+vLCVB4LKJwWN7mVh11ZRBHK1vsmFlu7hQEkeSqm5ufqKqWwBrTG8iyu9Xnl+5hzPys8jpndp+ARMTCof0YVhmGs8s393+xqbTCiVxrBCRR0XkHHd6BOeCuTER89H2A+w5fJwrJtnZRmciIlxZOIhlJYfYXmmNL7uqUBLHPKAIuBWni/UN7jJjIuZvS3fRu1sSF55qPeF2NldMysWXIDyzvNTrUEyEtDkeh4gkACtUdQxwb3RCMvHuQE0drxft49ppeaQm2b0bnU2/XqmcOzKL51aWcsfMEST6QuoSz3Qibf5F3XG914jI4CjFYwz/WFFKQ5PypdMGtb+xiUlXFQ6isrqOdzZXeh2KiYBQRgDsj3Pn+FLgaPNCVb0kYlGZuKWqPLV0F1PzMjilX0+vwzEn6NxR/cjskcLTy3ZzQYH1aNzVhJI4QhomtiUiMgu4H2fM8T+q6t1B68VdPwdnzPGvqupKd10JUA00AY2qWuguzwD+DuTh9KN1laoeOtEYTWxZsv0AJQeO8e0L8r0OxZyEJF8CVxUO5MF3t7Hn8HFybYz4LqWtvqpSReQ24EpgFLBYVd9tntrbsYj4cG4anA0UANeISEHQZrOBfHe6CfhD0PpzVXVCc9Jw3Qm8qar5wJvuc9NFPLV0N727JTF7TH+vQzEn6UunOTXcf/1op8eRmI7W1jWOJ4BCYB3OF/w9Ye57KlCsqtvdgZ+eBuYGbTMXeFIdHwHpItLeN8ZcN7bmGC8NMy4Tow7U1PH6+n1cPinXLop3AQP7dOe8Udn8fdlu6hqbvA7HdKC2EkeBqn5FVR8CvgCcGea+c4HAu4BK3WWhbqPAQhFZISI3BWyTraplAO6j9bXdRTy9bDf1TX6+NNXaYnQV150+hANH63lt3T6vQzEdqK3E8UlHhqp6Ih3PtNS5UHB/y21tM11VJ+Gc7dwcbv9YInKTiCwXkeWVldayI9Y1NPn585KdnHFKJvnZdlG8qzjjlEyGZqbx5JISr0MxHaitxDFeRKrcqRoY1zwvIlUh7LsUCGxPORAIHpS41W1UtfmxAqdb96nuNuXN1VnuY0VLL66qD6tqoaoWZmVlhRCu8dLrRfvYV1XLDdPzvA7FdKCEBOHLpw1m5a7DrN9zxOtwTAdpNXGoqk9Ve7lTT1VNDJjvFcK+lwH5IjLU7Zb9amB+0DbzgevEMQ04oqplIpImIj0BRCQNmIkzgFRzmevd+euBl0J+tyZmPb64hCF9u3PuSKt57GqunDyI1KQE/mIXybuMiN3S6VZv3QK8DmwEnlHVIhGZJyLNXZYsALYDxcAjwDfd5dnAByKyBlgKvKqq/3TX3Q3MEJGtwAz3uenE1pYeZsXOQ1x/eh4JCdZ9elfTu3sSl00cyAur9rC/ps7rcEwHCOU+jhOmqgtwkkPgsgcD5hW4uYVy24HxrezzAHB+x0ZqvPSnxSWkJfv4QqGN8tdV/duZQ3lq6S7+vGQnt88Y0X4BE9OsExnjqYqqWl5eu5crCwfRK9V66++qhmf14ILR2fz5o50cr7emuZ2dJQ7jqccWl9DkV776uTyvQzERdtNZwzh4tJ7nVlqvuZ2dJQ7jmSPHG/jLRzuZM7Y/eZlpXodjImxKXh/GD0rn0Q920OQPbplvOhNLHMYzf/loJzV1jcw7e7jXoZgoEBFuOnMYO/Yf5Y2N5V6HY06CJQ7jidqGJh5fvIOzR2QxJre31+GYKLnw1GwGZXTj9+9sw2kbYzojSxzGE8+uKGV/TT3fOMfONuJJoi+Bb5x9Cmt2H+bdLdajQ2dlicNEXX2jn4fe3cbEwemcNjTD63BMlH1h8kBy07tx/5tb7ayjk7LEYaLu2RW7KT10nFvPy8cZksXEk+TEBL5xznBW7TrM+1v3ex2OOQGWOExU1TY08bu3ipk0OJ1zRlofYvHqysKB9O+damcdnZQlDhNVTy/dRdmRWr47c6SdbcSxlEQf3zxnOCt2HmJx8QGvwzFhssRhouZ4fRMPvLON04Zm8Lnhfb0Ox3jsqimDGNA7lV++vgm/3dfRqVjiMFHz5JISKqvr7GzDAM5Zx3dmjmRt6RFeXVfmdTgmDJY4TFQcPFrP794u5pyRWUy1llTGddnEXEbl9ORXr2+mvtHvdTgmRJY4TFTc/8YWjtU38YM5o70OxcQQX4Jw5+xR7Dp4jL9+bON1dBaWOEzEFVfU8JePd3HN1EE2LKz5jLNHZDH9lL789q1iqmob2i9gPGeJw0TczxdspHuSj9susHEYzGeJCN+fPZpDx+q5b9EWr8MxIbDEYSLqnc0VvLmpgm+eewqZPVK8DsfEqDG5vfnyaYN5cslONpZVeR2OaYclDhMxx+ub+NFL6xmWlcbXzsjzOhwT4+6YOZJeqYn8+KX1dlNgjLPEYSLmd29vZffB4/zs0rGkJPq8DsfEuPTuydw5exTLSg7xwqo9Xodj2hDRxCEis0Rks4gUi8idLawXEfmNu36tiExylw8SkbdFZKOIFInItwPK3CUie0RktTvNieR7MCdma3k1D7+3ncsn5XK63exnQnTl5EFMHJzOz17dyIGaOq/DMa2IWOIQER/wADAbKACuEZGCoM1mA/nudBPwB3d5I/BdVR0NTANuDip7n6pOcKcFkXoP5sQ0+ZU7n19HWkqiNb81YUlIEO6+fBxVtQ38eH6R1+GYVkTyjGMqUKyq21W1HngamBu0zVzgSXV8BKSLSH9VLVPVlQCqWg1sBHIjGKvpQI+8v50VOw/x488X0NcuiJswjczpybfPz+fVtWUssDvKY1IkE0cusDvgeSmf/fJvdxsRyQMmAh8HLL7Frdp6TET6tPTiInKTiCwXkeWVlTZgTLRs2lfFvQu3MOvUHC6baLnenJh5Zw9nbG5vfvTiequyikGRTBwtdUYU3FSizW1EpAfwHHCbqja30fsDMByYAJQB97T04qr6sKoWqmphVpZ13x0N9Y1+bv/7Gnp1S+Rnl42x/qjMCUv0JfDrK8dTXdvIfz631lpZxZhIJo5SYFDA84HA3lC3EZEknKTxV1V9vnkDVS1X1SZV9QOP4FSJmRjw89c2srGsip9fPs6qqMxJG5nTk+/PGcUbGyt4bHGJ1+GYAJFMHMuAfBEZKiLJwNXA/KBt5gPXua2rpgFHVLVMnJ+qjwIbVfXewAIi0j/g6WXA+si9BROqBevKeHxxCTdMz2NGQbbX4Zgu4qufc46nu1/byJrdh70Ox7giljhUtRG4BXgd5+L2M6paJCLzRGSeu9kCYDtQjHP28E13+XTgWuC8Fprd/lJE1onIWuBc4PZIvQcTmu2VNXzvH2uZODid78+2VlSm44gIv/rCOPr1TOWWp1Zy6Gi91yEZQOKh7rCwsFCXL1/udRhdUnVtA1/4wxIqqmt55dYzyU3v5nVIpgtauesQVz/0EZOGpPPnG08jyWf3LkeDiKxQ1cLg5fbpmxPW2OTn5r+tYltlDb/70iRLGiZiJg3uw91XjOWj7Qf5yfwiu1jusUSvAzCdk6ryk/lFvLelkl9cMZbpp2R6HZLp4i6fNJAt5TU8+O42hmf14MYzhnodUtyyxGFOyG/fKuavH+9i3tnD+eKUwV6HY+LE9y4cScn+o/z0lQ2kd0viiskDvQ4pLllVlQnbI+9t595FW7h8Yi7fu3Ck1+GYOJKQINx/zQTOOCWT//jHGv65fp/XIcUlSxwmLE8uKeFnCzZy0dj+/PIL40hIsJv8THSlJPp46NrJjB+Uzq1PrWJhkSWPaLPEYUKiqjzwdjE/fqmIC0b3474vTiDRWrYYj6SlJPKnr06lYEAvvvHXlbxo3bBHlf3nm3b5/cpPX9nIr17fzKUTBvCHr0wmOdEOHeOt3t2T+Mu/ncbUvAxuf2Y1T3xY4nVIccP++02bauoamfeXFTy2eAc3TM/j3qsmWBt6EzN6pCTy+A1TOH9UNj+ZX8SPXlxPQ5Pf67C6PGtVZVpVsv8o//7kcrbvP8qPPl/A16bnWceFJuakJjnXPH75z0089N52iitq+N2XJlp/aRFkPx3NZ6gqz68s5eLffsD+mjqe/NpUbjxjqCUNE7N8CcL354zmnivHs2LXIWbd/z7vb7XhFCLFEof5lMPH6rnlqVV855k1jMzpyfxbzrCb+0ynccXkgbx083TSuyVx7aNL+Z9XNnCsvtHrsLocq6oygHMB/NkVu/nFPzdTdbyB/7hwJPPOHo7PmtuaTmZ0/168/K0z+NmrG/njBzt4bf0+fnrpqZw3ynpt7ijWyaFhWclB/nfBRlbtOsyUvD7899wxjO7fy+uwjDlpS3cc5AcvrGNrRQ3nj+rHf8wayagcO7ZD1Vonh5Y44tiqXYe4d9EW3t+6n6yeKdw5axSXT8q1axmmS6lv9PPY4h088HYxNXWNzB0/gFvOy+eUfj28Di3mWeKwxAFAQ5Off67fxxMflrB85yEy0pKZd/Ywrp2WR7dkn9fhGRMxR4418OB723h88Q5qG/ycOzKLG88YxvRT+tqPpVZY4ojjxKGqFO2tYv6avby0eg/lVXUMzujOdacP4eqpg+mRYpe6TPw4UFPHXz/exZNLSthfU8+Qvt25dEIul0/KZUjfNK/DiymWOOIscdQ1NrGi5BDvbq1kUVE52/cfJcknnJWfxZenDeacEf2snykT1+oam3h1bRnPrSzlw20HUIUxub04b1Q2543qx7jc3nH/P2KJo4snjiPHG1hbepjVuw6zctchPt5xkGP1TSQmCFOHZnDx+AHMHpNDevdkr0M1JuaUHTnOS6v38saGclbuOoRfISMtmUmD+1CY14fCIX0Yk9ub1KT4qs71JHGIyCzgfsAH/FFV7w5aL+76OcAx4KuqurKtsiKSAfwdyANKgKtU9VBbcXSVxOH3KweO1lN66BjbKo9SXFHDtsoaiitq2LH/6CfbDc9K43PDMzlrRBanD+9rVVHGhOHQ0Xre3VLJ+1v3s2LnQUoOHAMgQSAvM41ROT0Zmd2LYVlpDOzTjUEZ3embltwlr5NEPXGIiA/YAswASoFlwDWquiFgmznAt3ASx2nA/ap6WltlReSXwEFVvVtE7gT6qOp/thVLrCWOJr9S29BEbUMTdY1+jtU3ceR4A1XHGzgSMB0+1kBFdS37jtRSdqSWiupaGpr+9fdK8gl5fdMYntWDMbm9mDCoD+MG9aZXapKH786YrqWyuo4VOw+xYe8RNu2rZnN5NbsOHiPwq7Nbko+BfbrRr1cKfdNSyEhLJrNHMn17pJDeLYm0lER6pCbSIyXRmU9JJC3ZF/M9TLeWOCL5U3QqUKyq290AngbmAhsCtpkLPKlO9vpIRNJFpD/O2URrZecC57jlnwDeAdpMHCfqmeW7eX/rfvx+pcmvNKmi2jxPu8v9qtQ3+alr8FPX2ERtg5/ahiYa/aEl67RkH/16pZLTK5WpQzPI6e3MD0jvxvCsNAZldLcOB42JsKyeKcwak8OsMTmfLDtW38iug8coPXic3YeOUXroOKWHjlFZXcfaQ4c5UFNPdV37d6wnCCT5EkhOTCDZfUwKekzxJZCQAAkiJIgg0jyP+/xf8wkJuM//tezfzxxGwYCOvXclkokjF9gd8LwU56yivW1y2ymbraplAKpaJiL9WnpxEbkJuAlg8OATG9p0z6HjrN9zhARx+sJp/sP5EoSEBMEXsDzRl0BK4meXJyUmkJKYQGqSj9REH6lJzvwny9znvbslfWrq1S3JkoIxMap7ciKjcnq1eTNhbUMTB4/Wc+R4A0frGqmua+SoO1XXNnKsvon6Rj8NTX7q3Mfm5/VNfuob1X1swq/Q5PfjV/Cr4lentaRfFb/fWaafrAuchy9Oaejw9x/JxNFShV/wT+3WtgmlbJtU9WHgYXCqqsIp2+z2GSO4fcaIEylqjIlzqUk+BqR3Y0B6N69D6XCR/ElbCgwKeD4Q2BviNm2VLXers3AfKzowZmOMMe2IZOJYBuSLyFARSQauBuYHbTMfuE4c04AjbjVUW2XnA9e789cDL0XwPRhjjAkSsaoqVW0UkVuA13Ga1D6mqkUiMs9d/yCwAKdFVTFOc9wb2irr7vpu4BkRuRHYBVwZqfdgjDHms+wGQGOMMS1qrTmuNdsxxhgTFkscxhhjwmKJwxhjTFgscRhjjAlLXFwcF5FKYGcbm2QC+6MUzsnoDHF2hhihc8TZGWKEzhFnZ4gRYi/OIaqaFbwwLhJHe0RkeUstB2JNZ4izM8QInSPOzhAjdI44O0OM0HnitKoqY4wxYbHEYYwxJiyWOBwPex1AiDpDnJ0hRugccXaGGKFzxNkZYoROEqdd4zDGGBMWO+MwxhgTFkscxhhjwtKlE4eIPCYiFSKyPmDZr0Rkk4isFZEXRCS9lbIlIrJORFaLSER7SGwlzrtEZI/7+qvd8dlbKjtLRDaLSLE7Bns0Y/x7QHwlIrK6lbJR+SxFZJCIvC0iG0WkSES+7S7PEJFFIrLVfezTSvlofZatxRkzx2YbMcbacdlanDFzbIpIqogsFZE1boz/z10eU8dlWNQdL7srTsBZwCRgfcCymUCiO/8L4BetlC0BMj2M8y7gjnbK+YBtwDAgGVgDFEQrxqD19wA/9vKzBPoDk9z5nsAWoAD4JXCnu/zOlv7mUf4sW4szZo7NNmKMteOyxThj6djEGdG0hzufBHwMTIu14zKcqUufcajqe8DBoGULVbV5FPmPcEYX9FRLcYZoKlCsqttVtR54GpjbocG52opRRAS4CngqEq8dKlUtU9WV7nw1sBFn/Pq5wBPuZk8Al7ZQPJqfZYtxxtKx2cZnGQrPP8vm9bFwbKqjxn2a5E5KjB2X4ejSiSMEXwNea2WdAgtFZIWI3BTFmALd4lZbPNbKaWwusDvgeSmh/3N3pDOBclXd2sr6qH+WIpIHTMT5dZetzsiSuI/9WijiyWcZFGegmDk2W4gxJo/LVj7LmDg2RcTnVpdVAItUNaaPy/bEbeIQkR8AjcBfW9lkuqpOAmYDN4vIWVELzvEHYDgwASjDOd0OJi0s86J99TW0/Ysuqp+liPQAngNuU9WqUIu1sCyin2VrccbSsdlCjDF5XLbxN4+JY1NVm1R1As5Z5FQRGRNi0Vj5H/+UuEwcInI98Hngy+pWJAZT1b3uYwXwAs4pY9Soarl7sPmBR1p5/VJgUMDzgcDeaMTXTEQSgcuBv7e2TTQ/SxFJwvkC+auqPu8uLheR/u76/ji/+oJF9bNsJc6YOjZbijEWj8s2PsuYOjbd1zkMvAPMIgaPy1DFXeIQkVnAfwKXqOqxVrZJE5GezfM4Fy3Xt7RtpDQfUK7LWnn9ZUC+iAwVkWTgamB+NOILcAGwSVVLW1oZzc/Src9+FNioqvcGrJoPXO/OXw+81ELxqH2WrcUZS8dmGzHG1HHZxt8cYuTYFJEscVvIiUi35riIseMyLF5fnY/khHOKWgY04GTuG4FinDrD1e70oLvtAGCBOz8Mp/XCGqAI+IEHcf4ZWAesxTlQ+gfH6T6fg9OSZFsk42wpRnf5n4B5Qdt68lkCZ+Ccxq8N+PvOAfoCbwJb3ccMjz/L1uKMmWOzjRhj7bhsMc5YOjaBccAqN8b1uC28Yu24DGeyLkeMMcaEJe6qqowxxpwcSxzGGGPCYonDGGNMWCxxGGOMCYslDmOMMWGxxGGMMSYsljiMMZ4RkdEi8qCI/ENEvuF1PCY0ljhMVIiIisifA54nikiliLziZVwdTUQ+dB/TReSbkdp/G+vzJGDMlIDlbcbjljvujkvR4j7c7VRE7gl4foeI3BW0zUMiMr3dNwOo6kZVnYfTg22hW76bG0e9iGSGsh8TXZY4TLQcBca4XS4AzAD2RDsIcUTsuFfVz7mz6UCHJ46A/Ycrnfbj2aZOR3xtqQMub+cL/TScbuFDIiKXAB/g3D2Nqh534/C8TybTMkscJppeAy5y5z/Va6mIvOh2bV3U3L2125fQq+KMnLZeRL7oLv+O+3y9iNwWsI88cUaCe8Tdz0L312vz8t8DK4FHxR0pzi33MxG5NTDQ4F/dzb+sW3uNgO2ax124Gxju/nL+VdC+v9f8eiJyn4i85c6fLyJ/cee/Is6ocavdX/C+oP0jIj8SZ8TARSLylIjc4a7ytRBfq/G0RUSGicgqEZniLmoEHgZub2X70TjdYwxyY/uj+3f6q4hcICKLxRnx7pPOBFV1vpsQvxxqXMZjXvd5YlN8TEANTp89/wBScfoUOgd4xV3f3E9PN5z+fPoCVwCPBOyjNzAZp6+kNKAHTh9DE931eThfbBPc588AX3GX+4FpAdutdOcTcPoA6hsUbx6fHpHxDpzR71p8jcD32VL5oH1PA551598HluIM7vMT4OvAaOBlIMnd5vfAdUH7L3Q/w244I99tdWNs6zNoMZ7geJvngZE4fSxNCPo79sIZOa938+cSsP47OGOJNMcx1v2MVwCP4XQTPhd40d3+HOA3wEPAzUExlRClUThtCm9KxJgoUdW14gy2cw2wIGj1rSJymTs/CMjHSRC/FpFf4CSY90XkDOAFVT0KICLP4wzWs8otu0NVV7vzK3C+wD4AdqrqR24cJSJyQEQmAtnAKlU9EMZbaek1wrECmCxOz6x1OGdBhe77uBU4HydBLhMRcJJDcJfbZwAvqepxABF5uZ34Pggzxiyc3lqvUNWiwBWqWiUiT7qxHg8qdyFwA84wpztUdZ0bXxHwpqqqiKxzY0JV38HpZtx0IpY4TLTNB36N80uzL4CInIPT1fTpqnpMRN4BUlV1i4hMxukd9OcishBob2CmuoD5JpwvXXCusQT6I/BVIAfnl3CwRj5dlZsawmuERFUbRKQE5wv2Q5xeU8/FGSBpI07ieEJVv9/Gbloa4KdD4nMdwempdzrOWV2w/8NJeI9/EpBIdyBdVfe6PxAC4/AHPPdj3z2dml3jMNH2GPDfzb9EXb2BQ27SGIVTlYOIDACOqepfcJLNJOA94FIR6S7OGAqX4VT3hOsFnMF0pgCvt7C+HOgnIn1FJAVncKVwVONUIbXmPZxqnvdw4p8HrFZVxblI/AUR6QcgIhkiMiSo/AfAxSKSKs7odxfRtvbiCVaPMwb2dSLypeCVqnoQpxrsxoDF5wJvh/EappOyxGGiSlVLVfX+oMX/BBJFZC3wU/7VImcssFScsZp/APyPqq7EGWdhKc7Y0n9U1VWESVXrcb7knlHVphbWNwD/7b7GKzgD74Sz/wPAYvfCcEsXo98H+gNLVLUcqHWXoaobgB/ijIW9Fljkbhu4/2U4Z29rgOeB5ThnCScaT0tljuIkzNtFZG4Lm9wDBLaumo3ztzRdnI3HYeKS2yR3JXClqm71Op4TISI9VLXGrSJ6D7jJTawnsq88nOtIoY6F3dI+VgKnuUn3pLnVeYWqur8j9mc6jp1xmLgjIgU4o+292VmThuth92xsJfDciSYNVxPQ293fCVHVSR2RNNwm1KtxWpr5T3Z/puPZGYcxxpiw2BmHMcaYsFjiMMYYExZLHMYYY8JiicMYY0xYLHEYY4wJiyUOY4wxYbHEYYwxJiyWOIwxxoTFEocxxpiw/H+91N5ly2rglgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(np.sort(r), compressive_resistance_distribution)\n",
    "plt.xlabel(\"Masonry unit weight [kN/m$^3$]\")\n",
    "#plt.xlabel(x_label)\n",
    "plt.ylabel(\"Probability density\")\n",
    "#plt.ylabel(y_label)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Convergence assessment (see Ballio and Guadagnini, 2004)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {},
   "outputs": [],
   "source": [
    "realisations_list = []\n",
    "mean_list = []\n",
    "var_list = []"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [],
   "source": [
    "# def random_choice(sample_vector, sample_size):\n",
    "    \n",
    "#     r_sample_size = np.random.choice(sample_vector, size=sample_size, replace=False, p=None)\n",
    "#     r_sample_size_mean = np.mean(r_sample_size)\n",
    "#     r_sample_size_std = np.std(r_sample_size)\n",
    "\n",
    "#     mean_list.append(r_sample_size_mean)\n",
    "#     std_list.append(r_sample_size_std)\n",
    "    \n",
    "#     return r_sample_size, mean_list, std_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {},
   "outputs": [],
   "source": [
    "# random_choice(r, simulations_number)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {},
   "outputs": [],
   "source": [
    "def random_choice(sample_vector, sample_size):\n",
    "    \n",
    "    for i in range(sample_size):\n",
    "        # print(i+1)\n",
    "        r_sample_size = np.random.choice(sample_vector, size=i+1, replace=False, p=None)\n",
    "        r_sample_size_mean = np.mean(r_sample_size)\n",
    "        r_sample_size_var = np.var(r_sample_size)\n",
    "\n",
    "        realisations_list.append(r_sample_size)\n",
    "        mean_list.append(r_sample_size_mean)\n",
    "        var_list.append(r_sample_size_var)\n",
    "    \n",
    "    #return realisations_list, mean_list, var_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {},
   "outputs": [],
   "source": [
    "random_choice(r, simulations_number) #it takes 4m 32s on average"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = np.arange(1,simulations_number+1, 1)\n",
    "\n",
    "#print(n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAELCAYAAADDZxFQAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABECElEQVR4nO2deZwU1bXHv2cWGHZUEGVRRBnQaADFLT4jJq48txg1Kmrcn09NgtG4xCQuWdwiiURjoqAvRsXELbhFJa7BLaJCI7KKC8iOCgPD7Of90V1tdXVVd/VSPdv5fj7zme5bVbfOnaV+fc6591xRVQzDMAzDS1lrG2AYhmG0TUwgDMMwDF9MIAzDMAxfTCAMwzAMX0wgDMMwDF9MIAzDMAxfIhMIERkiIi+JyHwRmSciP0q0/1JEYiIyW0SeF5GBAdcfISILRWSJiFwZlZ2GYRiGPxLVOggR2R7YXlXfFZFewDvAccByVd2YOOeHwG6qeoHn2nJgEXAosBx4GzhFVT/IdM9+/frp0KFDiz0UwzCMDss777yzTlX7+x2riOqmqroSWJl4XSMi84FBnod8D8BPofYBlqjqUgAReQg4FsgoEEOHDmXWrFnFMN8wDKNTICKfBB2LTCA8BgwFxgBvJd7/GjgD2AAc7HPJIGCZ6/1yYN9orTQMwzDcRJ6kFpGewKPARCe0pKpXq+oQ4AHgYr/LfNp8Y2Eicr6IzBKRWWvXri2W2YZhGJ2eSAVCRCqJi8MDqvqYzykPAt/1aV8ODHG9Hwys8LuHqt6lqmNVdWz//r5hNMMwDCMPopzFJMBUYL6qTnK1D3eddgywwOfyt4HhIrKTiHQBTgaeiMpWwzAMI50ocxAHAKcDc0VkdqLtp8A5IjICaAE+AS4ASEx3naKq41W1SUQuBp4DyoF7VHVehLYahmEYHqKcxTQT/1zCMwHnrwDGu94/E3SuYRiGET22krodUV9fz7vvvtvaZhiG0UkwgWhHPPjgg+y1115ceumlNDc3t7Y5hmF0cEwg2hHr1q0DYNKkSXz3u99l8+bNrWyRYRgdGROIdoQjCL///e958sknOeigg1i5cmUrW2UYRkfFBKIdUVtbS1VVFT/60Y+YPn06CxYsYN999yUWi7W2aYZhdEBMINoRmzdvpkePHgAcddRR/Pvf/6a5uZnDDz+c+vr6VrbOMIyOhglEO8ItEABjxoxh6tSprFq1iunTp7eiZYZhdERMINoRXoEAOPTQQ9lxxx256667WskqwzA6KiYQ7Qg/gSgvL+fcc8/lhRdeYMmSJRmv//zzz6mtrY3SRMMwOhAmEO2IzZs3071797T2s846i/LycqZMmZLx+kMPPZSf/exnUZlnGEYHwwSiHeHnQQAMGjSIo446invvvZeGhobA61etWsWaNWuiNNEwjA6ECUQ7ora21lcgAM4//3zWrFnDE08EF71taGigqakpr3u/++67PPfcc3ldWyjTpk1j1KhRrXJvw+jMmEC0I4I8CIDDDz+cIUOGZExWFyIQN954I5dcckle1xbKqaeeSiwWI6r90w3D8McEoh2RSSCcZPWMGTNYunSp7zmNjY15C0RtbW2rr7UwgTCM0mIC0Y7IJBAAZ599NmVlZUydOtX3eENDQ95F/urr62lsbMzr2mLR0tLSqvc3jM6GCUQ7obm5mbq6uowCMXjwYP77v/+be++91/f65ubmvD2I+vr6vK8tFiYQhlFaTCDaCc76Bb9prm722msvVq5cmfYwdT79FyIQre1BWIjJMEqLCUQ7wankmsmDAOjatStA2nTXQgWirq6u1QXCPAjDKC0mEO0Ex4PIJhBdunQBSEsoO4JRSA6itUNM5kEYRmkxgWgn5OpBBAlEew4xmQdhGKXFBKKdUKhAdIQchAmEYZQWE4h2Qmt7EHV1dahqqz6kLcRkGKXFBKKd4AhEtllMTg7Cm6QuRogJKNiLaG5u5s0338zrWvMgDKO0mEC0E4oVYiokSe3uJ19uuOEG9t9/f15//fWcr3V7EJ988klBdhiGkR0TiHZC2FlMUYSYmpqaksJSqEDMmzcPyO8B73gQDz30EEOHDuXFF18syBbDMDJjAtFOiDoHoarMnDnTN87v7qvQqa5BIbAwOALheB9z584tyBbDMDJjAtFOCCsQ+eYgXn75ZQ488EDefffdtGNugcjFg7jlllt48MEHU9oqKyuT9ogIV111Vej+HPFyhKKszP58DSNK7D+snbB582bKysqSHkIQ+eYgPvroIwA2btyYdswrEC0tLaE8gHvuuYeHH344pc0RMMeeG2+8MWs/Do4wON/Ly8tDX2sYRu6YQLQTnEquIpLxvHxDTJ999hng7yF4Q0y33XYbe+yxR1abGxsb0wTJ8SDyyWWYB2EYpSWy/zARGSIiL4nIfBGZJyI/SrTfIiILRCQmIo+LSN+A6z8WkbkiMltEZkVlZ3shaD9qL/kKxIoVKwKP19XVJV83Njby8ccf8/HHH2e1xW//CXeIKVccYXBExwTCMKIlyv+wJuBSVd0V2A+4SER2A2YAu6vq14FFQKYg9MGqOlpVx0ZoZ7sg214QDvkW63MEIpsH0djYmPzKtnDNb/+JYiSp/UJMzc3NXHvttXz55Zc592sYhj+RCYSqrlTVdxOva4D5wCBVfV5VnafUm8DgqGzoSGTaj9pNvsX6wgpEU1MTDQ0NqGrWNRWZPIh8dqfLFGJ6/PHHue6667jsssty7tcwDH9K4qOLyFBgDPCW59DZwD8DLlPgeRF5R0TOz9D3+SIyS0RmrV27tij2tkVy9SCKGWLy8yCc15nwE4ggAQtDphCTEwbbsmVLzv0ahuFP5AIhIj2BR4GJqrrR1X418TDUAwGXHqCqewJHEg9PfdPvJFW9S1XHqurY/v37F9n6tkOhApEpxNTc3MyqVatSznPjzUGEFQi/EJOTZHf3GRavB+EOMTnHLC9hGMUj0v8mEakkLg4PqOpjrvbvA0cBEzQgkK2qKxLf1wCPA/tEaWtbJ0oPYs2aNcmHbtgQk7vPIPw8CEcwCvEg/EJMTlu2WV6GYYQnyllMAkwF5qvqJFf7EcAVwDGqWhtwbQ8R6eW8Bg4D3o/K1vZAWIGoqKgAghfKQXrROye8BMULMTk5Cq8H4bzPx4Pwhpi89wPzIAyjmET533QAcDrwrcRU1dkiMh64HegFzEi0/QlARAaKyDOJawcAM0VkDvAf4GlVfTZCW9s8Yae5ighdu3YN9CAgXQScNRAQbhZTGA8iKKTlPNyd2lK54A0xuYXC1kYYRvGpiKpjVZ0J+Pn7z/i0OSGl8YnXS4FRUdnWHgnrQQC+AuF+8Dc1NSWTxZDqQWTLQTQ1NYXyIIJCWs77TZs2ZR2HF2+IyU8gLMRkGMUjlECIyCBgR/f5qvpqVEYZ6YSd5gr+ApHJg8g1xOT0lUkggkp7OO+d2lJBLFu2DFVlhx12SLY5HoTTh7vvsCGmBx98kD322CPUSnDD6OxkFQgRuQn4HvAB4PxHKmACUSIaGhpoamoKLRBdunTJmIPwPrRXrFjB1ltvzeeffx56oZy3Ty9BISbnvVsgmpub0+oqOcLgnsMQxoPIJhATJkxI69cwDH/CeBDHASNUNfdpJ0ZRCFvJ1SFMiMnNihUr2HHHHUMJRK4hpjAexJYtW+jZs2fGMUFmgXAe+KUKMT311FOMGzculN2G0V4Jk9FbClRGbYgRTDEEIluIafDgwYhIqFpMxU5Sexe3BX26zxRiKmWSev78+Rx99NGcf37g+k3D6BCE8SBqgdki8gKQfOqo6g8js8pIIex+1A75CMR+++1HZWVlTiGmfHIQfiEm74wmv5Lj0HaS1I59H374YeT3MozWJIxAPJH4MlqJXD0IvxxEUIipoaGBtWvXMnDgwFACEXahXDYPwhticrNy5UrfPjNNcw1KUq9evZpHH32UCy+80PIOhpEjWQVCVf9SCkOMYMLuR+2QzYNwP1idh/GgQYOoqKgInMUkIqhqaA8i2zRXd9jK60G4Z1W5aWlpQVX517/+lTaOoBDTSSedxKuvvsohhxzC0KFDA+01DCOdrAFbERkuIo+IyAcistT5KoVxRpwocxDOwziTB1FXV5e8d6EhJr9V0F4PYs2aNb59trS0sGjRopS+GhoamDNnTmCSet26dUl73OMWEX784x8H2m8YRrgk9b3AncQL6x0M3Af8NUqjjFRaWyDq6+uTs3WKFWJy4/UgNmzY4NunqqYssGtubuaHP/who0ePZunS+GeWTElq79h+97vfBZ5rGEY4geimqi8AoqqfqOq1wLeiNctwE+U0V7dAZAoxOQKRa4gpKEkdZBsEC0RLS0vKuJqbm3njjTcAcEq9ewv4ffDBBxnvHUR9fT2LFy9OaTvyyCM555xzQvdhGO2dMAJRJyJlwGIRuVhEvgNsG7FdhotiJKmDchArVqygsrKSfv36ZfQg3CGmYnsQ3uKBboHwLpTzCoT3fu4Qk9dDyGUf7DPPPJPq6uqkLY8++ijPPvss99xzT+g+DKO9E0YgJgLdgR8CewGnAd+P0CbDQ5TTXDdt2kTv3r0RkYw5iKqqKioqKnLOQQQlqd14RcMtEO57qGpKcru5uTkpIH79uj0AEQktEMuXL+eZZ+Ilw/r27cuWLVs44YQTQl3r5dNPP+Wyyy5LEcGamhpuv/12m1VltHnCzGJ6G0BEVFXPit4kw0uxQkzl5eU0NzenPEwbGxuTJcIzhZi6du1KZWUl9fX1yYed14NoampK9uU8jFWVlpaWZOjHz4PIJBDucWTyIBy7g7ZBDVoE6KWuro4hQ4aktNXU1Piem+kBv2bNGlatWsURRxzBypUrGTNmTLLMx49//GOmTJnCLrvswhFHHJHVJsNoLcLMYtpfRD4gvqc0IjJKRP4YuWVGktraWrp06ZJ8+GYjyIPo1q0bkPpp2/1QzxRi6tq1KxUVFSkzjtznLliwgO7duzN//vzk/RzcD+0wIaYvv/wy5d4OqhooEE4uIkggWlpaQnkQ7733Xlqbd2ZUmMV4AwYMYNSoUclpxKeddlry2Pr164HsBQvzZe3atVx++eVZ9ww3jGyECTH9HjgcWA+gqnMA3+0/jWjIpdQ3BOcgHIHwfvKurIxXUskmEJWVlSkPNfe5sViMxsZGPv7447Rjfp/03WTyILwbHQUJxOeff+7bl/vcMB7Eueeem9YWJhQ0bdq0lCm4fixYsIBRo0Ylx5dPiOnpp5+mZ8+eKbO5pkyZwmOPPcbMmTN56aWXuOiii7jllluSYTLDyJdQH0lVdZnnU5N9NCkhuQpEUIjJyWFkCjFlykFUVlamTEl1P7ydTYecHEHQrKlsHsRTTz2Vsg7CPY5Vq1bx8ssvZ+zLb3W10x7mgex3jtfDcbd/9NFH7LTTTpx66qmUl5dnFKFf/vKXxGKx5PvXXnuNp59+mmXLliUX/2Xj6quvZvPmzSxevJgxY8YAcN5556Wcc8wxxwDBYmkYYQkjEMtE5BuAikgX4snq+dGaZbjJRyCc7T6dMtoNDQ2+AuENMWXKQVRUVKQIhFsEvAIRFGLK5EG8+eabHH300cBXXpBbIE499VTf69wEPcz33HNPjj32WN9jzz//PDvuuCMjRozwPe61+dNPPwXgnXfeYdiwYSxcuDDQHjfe0NTvf//7lPcrV65k4MCBPPjgg5xyyik59eV3zJLgRqGECTFdAFwEDAKWA6MT740SkY9AQOqn72LkIMJ4EE6OIhcPwmn74osvkm39+vVLG4ObqqoqX68g00N6+vTpvu2HH344I0eOZNq0ab7H77777pT3J554Ysr7ZcuWpbyfO3duoA2ZcNZsTJkyJa/rDaPYZBUIVV2nqhNUdYCqbquqp6nq+lIYZ8QJux+1g7OdqPsBHiYHERRiqq+v9w0xZfIg8gkxuRe59e7dO20Mbrp37x5qRlQ23A/zyy+/3Pec66+/PmMfXpH6+te/npMNQEq5EC/r1q1L5na893TvJ+7gJMHnzJnD+vXrOe6445I5GsPIhTCzmHYSkUki8piIPOF8lcI4I04xPIgwOYigEFNdXZ1viClTDiKfEJOfQGTzIIL6gnAhFvfDfOPGjXmFZQ499NDk64cffjjwPPfsLC+jR49Ovn7xxRdTjvXv35+ddtop7Zr77ruPwYMHp7XPnDkTgBtuuIHJkyczffp0/vCHPwTeOxOffvopV111lYWrOilhchD/AKYCTwL+AV4jUmpra9Pm5mciyhCTe/aMe62DU7IjHw/CaXPH1fv06ZM2BoeKiorkmo6gvvJh06ZNBT8ITzrppMBj2WZRZbu3++dz3XXXBYbMHBoaGtI8jZqaGjZt2sT2228PxMXowAMPZM899+Sss85KK2C44447AvGw2p577pnxfkbHI1SpDVWdrKovqeorzlfklhlJipWDCEpSZwoxOesH/HIQzrlffPFFUhj8BCKbB+GXWHY8CG/8H+JClk0gampqeP3119OOZ8JbLbbYPPfccxmPX3fddcnXb7zxBldffXXKynE32cTB4Ze//CUQ/zm+99579O7dm4EDByIiiAjf/va3+da3vsX777/PpZdemrzuvvvuS9mXIyj5D7BkyRKbMdVBCeNB3CYi1wDPk7qj3LuRWWWkUKhANDc309LS4puDaGxsTOYs/EJMTqjIKbXhF2Jyx8ELSVK7z3ME4m9/+1va+eXl5VkFwrk+V1ozlPLaa68lX3/jG98A8h+HH5dccolvuxOSchg3bhyvvPIKI0eOTLb99a9/ZdOmTYwbNy7Ztm7dOl5//XWOPfZYfvrTn3L11Vfzu9/9jiuuuCL0ok6jbRPmt7gHcDrxCq7OxwjFKrqWjHwWysFXD3DnYR0UYnI8C78Qk/MJNtNCObdA+OUgwghETU0NM2bMSLb16tUrcHzl5eW+az2CSoW0Z6ZOnVq0vl55JZzj75y3YMGCZNvkyZOZPHkyDz/8cLIu1W677Zasovvyyy/zm9/8BoB58+Zx4okn8p3vfAcgZa8OVS3JtrBGcQgTYvoOMExVD1LVgxNfJg4loqWlhdra2pxmMXk9COdhnU+IyenDEQj3A97rQYhIqBCT9wHR0tLCGWecwaRJk5JtmcYrInTv3p0tW7akfOLv1q1bThVb2wPekuNRc99992U8fuKJJyIifO9730uKA5ASzps2bRrHH388t956KxMmTKCsrIxRo0Zx//33U1ZWllw34jBnzpycxtnS0hJYEt4oLmEEYg7QN2I7jACckE0hISbnoZnPLCanDyfE5Mbpd9WqVUB8T4kwSWrHw3G3uVcYQ+bxOgLh3Wioe/fugTF7Ixzf/364Qs1///vfs55z2WWX8eCDDwLx6cSnn346ACeccAL77rsvIsKvfvUrRo8eTXV1NRs3buSxxx7j3nvvZcGCBey3337ce++9QFwUJk6cyD777MPgwYPp27cvX375JatWrUrxdIziEibENABYICJvk5qDOCYyq4wkue5HDbl7EJlmMXk9CAe3t7Fhwwaqqqro27dvqGmuXbp0SavS6qz4dnCP95JLLknZ28ERCO/c/h49egROizXaDu+//37y9c9//vPka2fmmpu33nqLHj16MGLECG677baUY9dccw2TJ08G4vuazJo1i1122YWbbrqJW265hf79+wfaUFNTQ69evaitraWhoYG+ffsWOKqOSRiBuCZyK4xAci31Dek5COd70EK5TLWYvDkIhx49eiT73bBhA71796aqqipUktrPg8gkEO77OnTv3p3ly5enhJhMIDom3/ve93zbHXGAuPfq5i9/+YvvNbNmzWLGjBlcddVV3H///ckqu/fffz/dunXjsMMOS+6e6LBq1SoWL17Mf/3Xf3W6/EmY/SDymtIqIkOI71+9HfHk9l2qepuI3AIcDTQAHwJnqeqXPtcfAdwGlANTVPXGfOwoJf/4xz+YO3duyqeiQslHIII8iKAktbuaa6YQk/tB3b179xQPok+fPlRVVWXNQQQJhDd85c5BeI+5Q0xVVVUp17jLkRuGl7FjxyZfu0uwu1+7cYoxOnz44Yf07NmTLVu2MG7cOE499VR+/etfR2dwKxMmB5EvTcClqrorsB9wkYjsBswAdlfVrwOLgKu8F4pIOXAHcCSwG3BK4to2zR/+8Ad+8YtfpMXTC6EYAhE2B5EtSe1+ULs9iI0bN2YUCK8gOfY5ZAsxeT0It0C4P9F1797dPAijqHhXsO+8884MGDCAoUOH8vHHH/Ob3/wmuaYk09fo0aP59re/zc4778z48eM59NBDefXVVznzzDN55pln2LhxY/IebWnVemSTlVV1JbAy8bpGROYDg1T1eddpbwJ+eznuAyxR1aUAIvIQcCzwgc+5bQJVZc6cOQDcfPPN3H///UXptxQehDsH4RTAcx68QTmIHj16pHgQvXv3plu3bslyEpmmufp5EO4yG97xBoWYamtr2WabbVLaTCCMtojzbABYunQpQLLEe1A4zI/x48dz22230a9fP3r37p32f1NsQvUuIt1ExL8WcrjrhwJjgLc8h84G/ulzySDAXSJzeaKtzbJy5UrWr1/PgAEDeOihh9KKq+VLrvtRw1c5CK9AdO3aFRHJuGGQ0+YQNgfheBB+OYjm5mZWr17Nli1bkklqN345iLAhJstBGJ2JZ555huHDh7PVVltRXl6e4qVEQZhifUcDs4FnE+9H51KsT0R6Ao8CE1V1o6v9auJhqAf8LvNp8/W7ROR8EZklIrPc87JLjRNWmjRpEmVlZdx6661F6bcQD8K7UM7ZtjRTiMl9PgRPc3V7EEEhJuePtqmpif3224+bb76Z5ubmooWYmpqaUmx1BMKmPRqdkaC90wshjAdxLfGQz5cAqjobGBqmcxGpJC4OD6jqY6727wNHARPUP+C2HHBXpxsMrPC7h6repapjVXVspmltUeO4kEceeSSnn346U6ZMSdkZLV+KOc3VTyC8ISbwFwivB+FNUjuzmNzTXN2zplavXs3q1atDz2LK5EG4j7uLBzpTX3fddddMPx7D6HC8+OKLGasP5EsYgWhS1ZyXLUr84+NUYL6qTnK1HwFcARyjqrUBl78NDE+UGu8CnAy06RLjsViMHXbYga222oqf/OQn1NfX511i2U0xcxCVlZWhBMJ93AkZeWcxOSGmlpYWampqfD0Id86joaGBhoaG0CEm93j9QkzbbbcdQEpBuVzCcIbRVtlmm20YPnw4P/vZz5g5cyYrVqygpaUlmR9saWlBVVO+Dj744EhsCZOkfl9ETgXKRWQ48S1Hw5TJPIB4Dae5IjI70fZTYDLQFZiRCEG8qaoXiMhA4tNZx6tqk4hcDDxHfJrrPao6L5eBlZo5c+Yk9xYYOXIkxx13HHfccQeXX355Qcqej0A4sUnvLKYuXbqkFblrbGxMKbXhPh++clt79eqV8qB2PAinRHafPn2oq6tLEQhnCmpDQwPNzc3U1dWhqmkC0dLSklOSWkSSZajdmEAYrcVuu+3G4MGDufTSSzn88MOBeJXbbt26MXDgQCZOnMixxx5Lc3MzPXv2ZN999wW+qk+lqlkTzk7ItpRrMcIIxA+Aq4mvop5G/KH9y2wXqepM/HMJzwScvwIY73r/TNC5bQ0n7n3ccccl26644goef/xx7r777rQa+7mwefNmRCT5aTwMIkLXrl3TFsp5Q0yqmjXE5BYI57hTLK+hoSFZE6d3795s2LAhWR/J7UE4ouF4I8XwIIYOHZo27q233jr7D8cwfBg9ejSzZ8/mzjvvZOnSpcyaNYuKigoeeOABGhoaeOedd/jmN79Jnz59kp/mvX/HDn5Rc+/+426iTDIXSpiFcrXEBeLq6M1pn8yfP5/m5uaU3cn23Xdfxo0bx6RJk7j44osD/5iy4Ww3musfkLvaaVCIyanxnynEtHHjRrp27UqXLl1SZjt16dKFxsbGpED06dOHtWvXJsXBnYNw8ijOd2+S2k8g3F6D3zRXJ8TkZtSoUVl+KkZH5bnnnuOQQw6hrKyMkSNHsnDhQh555BH2228/Bg8ezDbbbMPWW2/NnDlzEBE+/PBDKisrqa6uDtX/oEFfTaJ0ys13BgJ9GhF50r3FqPerlEa2dZwEtXcv4iuvvJLPPvuMBx7wm6gVjlz3o3bwEwivB+F8zxZickJkznFHLBoaGpILfJwcBMQ9BrcH4XgOQR6EN8Q0ffr0FEH08yDKysrS/klHjMh7JrbRBnHKhV955ZXJtmOPPZYlS5YQi8WSf2dbtmzhsMMOS/4NzZo1i2nTpvHd736XQYMGoaqsW7eORYsW0a1bN6qqqvja174WWhw6M5k8iN+WzIp2TiwWo6qqiuHDh6e0H3bYYYwePZqbbrqJ73//+3ktasl1LwgHt0AE5SCc9mwhJkcgMnkQziwmiAtEbW0tO+ywA/CVMDgehF+Iyb1j2THHpNaB9MtBOP24S2u4F80ZbYvRo0fzj3/8gzVr1rDPPvukHFNVbr/9dn7wgx/w6quvMnLkyGShvbq6Oqqqqth999057bTTOOOMM9h5551Trvd+gOjZsycnn3xytAPqJAQ+sTzbi74BfAF8DrxhW46mMmfOHHbfffe0T7QiwhVXXMHChQt54on8nK7a2tq8BKJLly5ZcxDO90whppqamuSuZm6BqKyspKWlhS+++AKIexBuj+HLL7+kX79+yffu734eRKZ9HPxCTH7tQecZueOuWRTEoYcemtY2fPhwRo8eDcApp5ySbH/99dfZcccd2XvvvZk2bRqfffYZv/3tb5ObE1100UXMnTuXAw88MKUKq/OhY8KECSxfvpzjjz++kGEZORJmodx/Ey+qNxm4HVgiIkdGbVh7wSmxERT/PuGEExg2bBg33HBDXjVWiuFBBOUgwoSYNm7cGBhigvi2k5AaYqqpqWHz5s1pAuHnQVRWVrJ48WI++uijwLH4hZjcdhtf8emnnwYec0o7ONxzzz3JrU0BzjzzzOTrN998M+1691TKc889l2effZYtW7bwt7/9jTvvvJOamhoWLlzIe++9xxdffMF9993HunXrWL16dcoki5NPPpmBAwdy6aWX8s1vfhOI/0533333jGNz5wGM0hAm5nErcLCqjlPVg4CDgd9luabTsGrVKtatW5eWf3CoqKjgsssu4z//+U/oLR/dFDvE5BaIQkNM8JVAuENMzgZCjkBkSlJXVlbyyiuv8MknnwSOJSjEtH79+rRzM4XxgvZkbk8ccsghydeLFi1KO55JNJ2HMcDbb7/NWWedlfI3eeeddyZfl5eXc8cddyTfqyovvvgin3zyCRs3buTuu++mrKyMqqoqTjrpJC644AJ69uyZ/N307duXiooKttlmG7bddtv8Bmu0OmEEYo2qLnG9XwoUvkS4g+CU2Mg0g+bMM89k22235aabbsq5/2J6EN4cRNgQk1cgvB5EWVkZPXv2DBSITCGmMF6A30rqIDLNLrn55ptTHnqtzbhx43I6f+HChSnJ+6233pobb/yqCv4pp5yS9vM88MADk68rKyv57W9/yzXXXJMMIVVUVLBs2TI++eQTqqqqUjyKCy+8kCeeeII33ngj2bbDDjtEsmLXaJuEEYh5IvKMiJyZKJHxJPC2iBwvIp0+IBg0g8lNt27dmDhxIs8++yyzZ8/Oqf98ZzG5d23LFmLKVIvJnYNwC4nbg+jduzcikiYQTtLYG2JyexBhHv5BHoQfmY6VlZVx4YUXpu1MFiWZKnVOmzYt5f2QIUMCzoxTXV2dFqZ0ewXO9p4OS5Ys4emnn2bXXXdN7hl96aWXcu2116acN3jw4OSEghkzZrB8+fLksaOPPpr99tsvo11GxyWMQFQBq4GDgHHAWmBr4pv+HBWZZe2EWCzGkCFD2GqrrTKe97//+7/06tWLm2++Oaf+C/Eg3MX6KioqEBHfEJO3mmtQDiLIg3AExIkzO+UvvB6EO9TlEMaDKJZAOMd++MMf+h53b2vqcNddd2W1LxPu2VlB9gCcffbZzJo1K/DcJ598EkidqbXNNtskQ2qOp+DeQGnnnXemV69efPDBB+y///6h7O3evbvF+o0kWQVCVc/K8HV2KYxsy7hLbGSib9++XHDBBfztb39L1oMPQ76zmLwhJuehnMssppaWFjZt2uSbg3ALhLOXcLYQk0OQQOy0006+D8lcQkz5igekr4CdPXt23gscITWmn8mefv36MXXq1OTv2W/V/FFHxT+LOSE0Z7+Rvffem5/85CfJtTb5/K0YRhBhZjHtJCKTROQxWyiXilNiI+wK3okTJ1JRUcFvfxt+iUmxchDuMhlhQ0xOHahMIaa1a9cmjxcqEAcddBB77bVX2lhy8SAKwS0Qe++9N6NGjUqOLYhMlWOrq6uZMGEC1113ne/5W2+9NYMHD+aPf/wj8FXo7Re/+EXynNNOOy3l78UrYmVlZdx8881Zw1OGkQ9hPpr9g3hV1ieJ7y1tJFiwYAFNTU2hPAiIb6x+xhlncO+993LNNdcwYMCAjOc7VVDzXQcR1oMICjG56zC5j7tDTKtWrUoLbzghJqc2kpN7cPDOYvJ77SZomqsfheyw5RUIIKW+VlhOOOEE1q1bx/77709lZSW/+MUvuOaaawDSVogvW7Ys5b1jw1VXxXfi/etf/+prY1ut3WN0LML8N9Wp6mRVfcmzeK7TEyZB7cUpBT558uSs5+ZTydXBm4PwE4igaa7OcaeMRqZprgDDhg0DUgWioqKC7t27U1FRkeZBuD+Vux/+QaGkYuUg3PjljA466CDf/nbbLXg7dL+1LT//+c956aWXfENFhW4R6cwy8q7aN4woCPPXepuIXCMi+4vIns5X5Ja1A4JKbGSiurqa448/njvuuCNlo3I/ChWIXHMQ3hCT14PwWygHXwmE8+D/8ssv6du3bzIp7hUI9zRJdz/u6bVuolgQN2vWrJR9w1WVsWPHcsEFF6Sd6xaBiy66KGvfmT4wFCoQF198MYsXL056OIYRJWH+WvcAzgNuJL5o7lasThMQ9yC+9rWv5ZREhXgp8A0bNjB16tSM5+WzH7VDPjmIoBBTUKkNB0cgunfvnhRL5xN6eXl5WojJ7UG4+wkqt5FLiCmsBzFs2DAmTJiQ1u5M6XQS7168OYlcV8cfeWRhRQhEhF122aWgPgwjLGEE4jvAMFU9SFUPTnx9K2rD2gOxWCyvEtN777031dXVzJw5M+N5xfIgvCEmb7E+bw7CEY6gHIQ3xLTTTjslX48ZMwaIz9py7pfJg3A//IMEwutB7LHHHkHDLjg2f9pppzFp0qSURLG3OJzDRx99lLNA/PrXvy7IPsMoJWEEYg7QN2I72h2rVq1izZo1OeUf3Dg16zORz37UDt5iffmEmLw5CL8QU1lZWXKRFcCee8ajj45AlJeXZxQI9wPdsdeL14PwJm7dFCoQ5eXlXHLJJSnrCe6//37uv/9+zjvvPK644opk+9ChQznppJNy7j8smYQwE5MnTy7pYkCj4xImNjIAWCAibxPfVQ4AVT0m+JKOT5gSG5morq7mueee890sx6FQD6KlpSU5E8pdkC/XEFMmD2KHHXZI+YTvFQjvHtju/iB1IVmYEFPv3r3p2bNn2jlOTN4rEOecc07WUF42+vTpw4QJE3xDUtdddx2/+tWvsvbx85//nD//+c8AXHbZZSkC5MdHH32Ud/nyH/zgB3ldZxhewgjENZFb0Q7JZwaTm+rqaurr61m2bJnv9plQuEBAfK2GO8TkzkGEDTF5cxBuD8LJPzj4hZi8uOP4YQSivLw86QWdd955acdffPHFZIlpRyCmTZvGunXruPjiiwsWiEy4k84zZ84MXFh3/fXXc/311wNwyy23ZO036G/CMEpJmC1HbUqrD7FYjMGDB+e9D7Kz+9nChQsjF4iGhobkJ1Z3DiJMiKm8vDzlWkhNUnsFol+/fpx++unJjdv9vCP39E+3QASFmMrLy6msrKS+vj7FW7nhhht46qmnUspQO3zta19Lhmj22WefrKWkc2HTpk2+7QcccEDR7mEYbYGsAiEi+wF/AHYFugDlwGZVzbzEtIMTtsRGEM52h4sWLUo+TL0UIhDOJ1lHINyroXMJMfXq1Stt/4UuXbok+/cKBMB9992XfO3nQbjDQI5Yue/rxREZ76fzK6+8MmU7Snff7k/2b731lm+/+WLlLIzOQpgk9e3AKcBioBtwbqKt09LQ0MD8+fPzzj8ADBgwgF69evnW9HcodJorxG3NloNwJ5zLyspSQkx+axYqKysZOHAg3/72tzniiCMy2pEtKRvWgwiLIxCdZVN5w4iSUBP4VXWJiJSrajNwr4i8HrFdbZpcS2z4ISKMGDEi40ymKHIQmVZSO6/dHoTfqufKykqqqqrSdijzI9saEbdAnHrqqb7n5LK4zM+DMAwjP8IIRK2IdAFmi8jNwEqgU/vYhSaoHaqrq3nttdcCj9fW1lJRUZFXRVFvDsKdpA7KQUD84e/OQfh5ELnYk20VtGPLPffck7LlpZtcpq6aQBhG8QjzX3R64ryLgc3AEOC7URrV1onFYnTt2jWZR8iXESNG8Omnn6atE3DIt5IrpOcgwkxzhfgDPUyIKSzObKYgHA+ie/fuRSlAZyEmwygeYWYxOZsF14nIZGCIZwvSTke+JTa8ODuELVmyxHdRVCEC4c5BZAsxuR/43hDTwIEDk8f69OnD//zP/3DYYYeFtiPbLC/HgyhWvaVcPYj33nsv46Y+htGZCTOL6WXgmMS5s4G1IvKKqv44WtPaLrFYjPHjxxfcj3smU1QC4Q0xhfEggnIQZWVl/OlPf8rJjmw77TkP59YSCGf9hGEY6YT5L+qjqhuB44F7VXUv4JBozWq7rF69mtWrVxecf4CvBCIoUZ3vftRQWA7CXe670A3qs3kQxRYIBwsxGUbhhBGIChHZHjgJeCpie9o8hZbYcNOzZ08GDhwYONW1WB5EY2NjSg5CVWlpackYYlLVtBxEPng9CO+DO6oQk2EYhRNGIK4HngOWqOrbIjKM+JqIjIjIEBF5SUTmi8g8EflRov3ExPsWERmb4fqPRWSuiMwWkeDd3EtMsWYwOYwYMSISgXA8BmcdhDvEBHHvwfEU3A9tJ8RUV1dHc3Nz0QXCm7eJKsSUa5VVwzDSCZOkfhh42PV+KeFmMTUBl6rquyLSC3hHRGYA7xMPV/05RB8Hq+q6EOeVjFgsxqBBg/IupOalurqahx9+2PdYbW0t22+/fV79Oh5EbW0tLS0tgQJRXl6e8qnbCTF56zDlizfE5BUC50HengVi8uTJgavADaM9U9g0nAyo6kriayZQ1RoRmQ8MUtUZ0H5DAYWW2PBSXV3N559/zvr169NEpxghJqdukHvDIPhKILwPZifE5C31nS/ZPIiOEGKy6qlGR6Ukq4lEZCgwBsilKI4Cz4vIOyJyfiSG5YhTYqOYAuEu2uelGALheAJeD6K5uZnGxsa0B7YTYvKW+s4Xrwfh3T8hU4jpJz/5Scb9oP0499xzgezrLwzDyE5WgRCRncK0Zbi+J/AoMDExGyosB6jqnsCRwEUi8s2A/s8XkVkiMmvt2rU5dJ87CxcupLGxsSgJagf3VFcvxchBBAmE40H4CURTU1PSgwjaejMsbg9i5cqV3HHHHSnHHYHwW1Ny8803M2/evJzud80111BfX++7Z4RhGLkRxoN41KftkTCdi0hl4voHVPWxXAxT1RWJ72uAx4F9As67S1XHqurY/v3753KLnCl2ghridf8rKirSBEJVizLN1RtiyiYQ3hBToTkIt0Bst912JQkx5VOaxDCMdAJzECIyEvga0EdEjncd6g1k3g4rfr0AU4H5qjopF6NEpAdQlshd9AAOIz6bqlWJxWJ06dIlGRYqBpWVley8885pIaa6ujpUtWg5CPc6CIgLhHv6q9ueLVu2sGHDBqBwgQhbaqPY6yAMwyicTEnqEcBRxPejPtrVXgOkb+uVzgHE6zjNFZHZibafAl2J7y/RH3haRGar6uEiMhCYoqrjiW9z+ngi4VgBPKiqz4YdVFQUq8SGl+rq6jQPopD9qCEuBOXl5RlzEEEhppqamqJ5ENl+VsX2IAzDKB6B/72qOh2YLiL7q+obuXasqjOBoCklj/ucvwIYn3i9FCheoL9IxGKxrPsf5MOIESN4/vnnaWlpSZaIKKTUt0OXLl3SPIhcQ0yF5iAchgwZ4ttuHoRhtF3CfBReIiI/BYa6z1fVs6Myqi2yZs0aVq1aVdT8g4OzP/Wnn36a3H60GALRtWvXjDmIoBBTY2MjGzZsoLKyMhmqKoTPPvsscBwmEIbRdgkjENOBfwP/ApqznNthKWaJDS/umUzFFghviMm7DiLTLKbevXsXZV2BuyKsFwsxGUbbJYxAdFfVKyK3pI0TxQwmB/daCKeUdlQCkS0H4Q4xFSu8lAnzIAyj7RJmmutTIlJ4bet2TiwWY+DAgfTr16/offvtT13IftQO7hxE2Gmu7hBToQnqMHTr1i3FPsMw2g6ZprnWEF/NLMBPRaQeaEy8V1WN/unRhih2iQ03IpI2k6lYHsTKlSsB/yS1Xw7C2S/CCTFFzRtvvMHTTz9t5bkNow2SaRZTYTUWOhCNjY188MEHHH744ZHdY8SIESn7Uxc6zRXiAuEITS45CCfENHjw4LzvnY3LL7+curo6dt99d3bffffI7mMYRv6E2VFuT5/mDcAnqtpUfJPaHgsWLCh6iQ0v1dXVTJs2jS1bttCtW7eieRAOuayDcEJMudZByoWbbropsr4NwygOYXIQfwTeBO5OfL0JPAQsEpHwmxO3Y5wZTFGFmCDuQagqH374IVC8EJND2GmupQ4xGYbRdgkjEB8DY1R1r8R2o6OJ7+lwCHBzdKa1HaIoseHFu/1osZLU3tdhk9SlmsVkGEbbJcw015GqmiypqaofiMgYVV3aXvd0yJU5c+aw2267RTrTZvjw4cBXVV03b95MVVVVcmV1PviFmMLkIOrr64HCy2wYhtG+CfP0WSgid4rIQYmvPxIPL3UlPqupwxOLxSLNP0B834WBAwemeBCFhJcg/xCTgwmEYXRuwgjEmcASYCJwCbA00dYIHByRXW2GtWvXsnLlykjzDw7uqa7FFohcktQOFmIyjM5NmD2ptwC3Jr68bCq6RW2MKEtseBkxYgSPPBLfaqO2trZggcg3B+FgHoRhdG4yLZT7u6qeJCJziS+YS0FVo/9I3QaIssSGl+rqatavX8/69esjCzFly0FYiMkwDIdMHsSPEt+PKoUhbZVYLMb2229P1LvVQWrRvtbKQViIyTAMh8AchKquTHz/JNE0PPF6DfB5CWxrE0RZYsOLM4222AJRWVmZrMqaSw7CPAjD6NxkTVKLyHnE96D+c6JpMPCPCG1qMzglNkolEM7+1AsXLixoP2oHJ+/gzkWE2TDIwQTCMDo3YWYxXUR8+9CNAKq6GNg2SqPaCgsXLqShoaEkCWqIf3ofNmxYJB6EQy4hJhMIw+jchBGIelVtcN6ISAU+SeuOSClKbHgZMWIECxcuLMosJkcg3B5EmIVyzrXF2E3OMIz2SxiBeCWx5Wg3ETkUeBh4Mlqz2gaxWIzKykpGjhxZsntWV1ezePFiNm3aFIlAhNkwCMx7MAwjnEBcCawF5gL/AzwD/CxKo9oKpSix4WXEiBHU19dTX19fkhBTkAdhAmEYRphaTOOAB1T17ohtaXPEYjEOOeSQkt7TmeoKhVVyhcxJ6oaGBlpaWgJzEDbF1TCMMAJxJvAnEVkP/DvxNVNVv4jSsNZm3bp1rFixoqT5ByiuQGTKQTgF+SzEZBhGEFlDTKp6hqpWA98FlgN3EA85dWhKWWLDzXbbbUevXvHN/Aqd5popBxEkEBZiMgzDIcw6iNNE5M/E10IcAtwOHBi1Ya1NKUtsuHH2p4bieRDuMFJZWRkiwpYtW9KOud9biMkwjDAhpt8DHwJ/Al5S1Y+jNKitEIvF2G677dh229Iv+aiuruadd96JJAcBca+hrq4u+dp7DMyDMAwjXIipH3A2UAX8WkT+IyJ/jdyyVqaUJTa8OCU3oshBQDwPESQQFmIyDMMhTIipN7ADsCMwFOgDtERrVuvS1NTEvHnzWk0gdt11VwD69u1bUD9+ISZI9SAsxGQYRhBhQkwzXV+3q+ryaE1qfUpdYsPL8ccfz/Tp09ljjz0K6ifIg7AQk2EYYQgTYvq6ql6oqg/mIg4iMkREXhKR+SIyT0R+lGg/MfG+RUTGZrj+CBFZKCJLROTKsPctBq1RYsNNRUUFxxxzTLICa77kk4PYaqutKCsrY9CgQQXd2zCM9k8YDyJfmoBLVfVdEekFvCMiM4D3geP5qjpsGiJSTnw67aHEp9a+LSJPqOoHEdqbpDVKbERBphyEM4vJKxCDBg3igw8+YPjw4aUx0jCMNktkApHYT8LZU6JGROYDg1R1BpDt0/E+wBJVXZo49yHgWKAkAjFnzhx23XXXtAdreyOfHAR8lSQ3DKNzEyZJfUCYtix9DAXGAG+FvGQQsMz1fnmiza/v80VklojMWru2OOv3YrFYq+UfikmmHETQQjnDMAyHMMX6/hCyzRcR6Qk8CkxU1Y1hL/Np8y0xrqp3qepYVR1bjG1B169fz2effdZq+Ydikk+S2jAMwyHw6SAi+wPfAPqLyI9dh3oD5WE6F5FK4uLwgKo+loNdy4EhrveDgRU5XJ83rVViIwocYfCGkdw5iFJWqjUMo32RyYPoAvQkLiK9XF8bgROydSzxJMNUYL6qTsrRrreB4SKyk4h0AU4Gnsixj7xorRIbUWAehGEYhRD4dFDVV4hvFvR/qvpJHn0fAJwOzBWR2Ym2nwJdiYeo+gNPi8hsVT1cRAYCU1R1vKo2icjFwHPEvZV7VHVeHjbkTCwWY8CAAQwYMKAUt4uUsrIyDjroIPbcc8+UdstBGIYRhkwhpt+r6kTgdhFJi/+r6jGZOlbVmfjnEgAe9zl/BTDe9f4Z4psTlZTWLLERBS+//HJaW0VFReA0V8MwDIdMTwen3tJvS2FIW8ApsXHxxRe3timRUl5eTkNDfJtxy0EYhhFEphDTO4nvr5TOnNZl0aJF1NfXd4gEdSbcXoN5EIZhBJH16ZBY83At8WJ9FcTDRqqqw6I1rfS0domNUmECYRhGGMI8HaYClwDvAM3RmtO6zJkzh4qKimQ11Y6KWxQsxGQYRhBhBGKDqv4zckvaALFYrEOU2MiGsy81mAdhGEYwYZ4OL4nILcBjQL3TqKrvRmZVKxGLxRg3blxrmxE5FmIyDCMMYZ4O+ya+u0tzK/Ct4pvTenz++ecsX768w+cfwEJMhmGEI6tAqOrBpTCktelIJTayYR6EYRhhyLRQ7seeJgXWATNV9aNIrWoFOlKJjWxYDsIwjDBkqsXUy/PVm3iY6Z8icnIJbCspsViM/v37d4gSG9kwD8IwjDBkWih3nV+7iGwN/At4KCqjWoM5c+YwatSogrf5bA9YDsIwjDCE2Q8iBVX9nOAaS+0Sp8RGZwgvgXkQhmGEI2eBEJFvAV9EYEursXjxYurq6jpFghpMIAzDCEemJPVc0ndx25r4xj1nRGlUqeksJTYcnCS1iKQkrA3DMNxk+vh4lOe9AutVdXOE9rQKnaXEhoPjNZj3YBhGJjIlqfPZJKhdEovFGDlyZHIHto6OCYRhGGHIOQfREYnFYp0m/wBfCYPNYDIMIxOdXiAaGxvZa6+9OkUNJgcn72AehGEYmej0T4jKykoefzxtB9QOjYWYDMMIQ6f3IDojJhCGYYTBBKITYjkIwzDCYALRCbEchGEYYTCB6IRYiMkwjDCYQHRCLMRkGEYYTCA6IeZBGIYRBhOITojlIAzDCIMJRCfEPAjDMMJgAtEJsRyEYRhhiEwgRGSIiLwkIvNFZJ6I/CjRvrWIzBCRxYnvWwVc/7GIzBWR2SIyKyo7OyPmQRiGEYYoPYgm4FJV3RXYD7hIRHYDrgReUNXhwAuJ90EcrKqjVXVshHZ2OiwHYRhGGCITCFVdqarvJl7XAPOBQcCxwF8Sp/0FOC4qGwx/LMRkGEYYSpKDEJGhwBjgLWCAqq6EuIgA2wZcpsDzIvKOiJxfCjs7CxZiMgwjDJE/IUSkJ/AoMFFVN4pI2EsPUNUVIrItMENEFqjqqz79nw+cD7DDDjsUy+wOjQmEYRhhiNSDEJFK4uLwgKo+lmheLSLbJ45vD6zxu1ZVVyS+rwEeB/YJOO8uVR2rqmP79+9f7CF0SJwchIWYDMPIRJSzmASYCsxX1UmuQ08A30+8/j4w3efaHiLSy3kNHAa8H5WtnQ3zIAzDCEOUHsQBwOnAtxJTVWeLyHjgRuBQEVkMHJp4j4gMFJFnEtcOAGaKyBzgP8DTqvpshLZ2KkwgDMMIQ2RPCFWdCQQlHL7tc/4KYHzi9VKg82wSXWJMIAzDCIOtpO6EWA7CMIwwmEB0QsyDMAwjDCYQnRATCMMwwmAC0QmxldSGYYTBBKITYrWYDMMIgwlEJ8RCTIZhhMEEohNiAmEYRhhMIDohloMwDCMMJhCdEPMgDMMIgwlEJ8SS1IZhhMEEohNiISbDMMJgAtEJ6devH0cffTT7779/a5tiGEYbxmIMnZDKykqeeOKJ1jbDMIw2jnkQhmEYhi8mEIZhGIYvJhCGYRiGLyYQhmEYhi8mEIZhGIYvJhCGYRiGLyYQhmEYhi8mEIZhGIYvoqqtbUPREJG1wJfABldznwzv3a/7AeuKYIb3fvmeG3TMrz3TGL3vbcyda8zFGm+QTfmcV6wxR/07DrIpn/Pa8ph3VNX+vkdUtUN9AXeFfe95PSuK++d7btAxv3Ybs405aMzFGm8uY852XrHGHPXvuLOO2f3VEUNMT+bw3nssivvne27QMb92G7ON2fu+Ncec7bxijTnq8ebSb0cac5IOFWIqBBGZpapjW9uOUmJj7vh0tvGCjbmYdEQPIl/uam0DWgEbc8ens40XbMxFwzwIwzAMwxfzIAzDMAxfTCAMwzAMX0wgDMMwDF9MIHwQkR4i8hcRuVtEJrS2PaVARIaJyFQReaS1bSkVInJc4nc8XUQOa217SoGI7CoifxKRR0Tkf1vbnlKR+J9+R0SOam1bSoGIjBORfyd+1+Py7afTCISI3CMia0TkfU/7ESKyUESWiMiViebjgUdU9TzgmJIbWyRyGbOqLlXVc1rH0uKR45j/kfgdnwl8rxXMLQo5jnm+ql4AnAS026mgOf4/A1wB/L20VhaXHMeswCagClie902jWH3XFr+AbwJ7Au+72sqBD4FhQBdgDrAbcBUwOnHOg61teynG7Dr+SGvb3QpjvhXYs7VtL9WYiX/oeR04tbVtL8WYgUOAk4l/EDiqtW0v0ZjLEscHAA/ke89O40Go6qvA557mfYAlGv/03AA8BBxLXHEHJ85ptz+jHMfcIchlzBLnJuCfqvpuqW0tFrn+nlX1CVX9BtBuw6c5jvlgYD/gVOA8EWmX/9O5jFlVWxLHvwC65nvPinwv7CAMApa53i8H9gUmA7eLyH9TguXsJcZ3zCKyDfBrYIyIXKWqN7SKddEQ9Hv+AfFPl31EZBdV/VNrGBcRQb/nccRDqF2BZ0pvVqT4jllVLwYQkTOBda6HZ0cg6Pd8PHA40Be4Pd/OO7tAiE+bqupm4KxSG1Migsa8Hrig1MaUiKAxTyb+YaAjEjTml4GXS2tKyfAdc/KF6v+VzpSSEfR7fgx4rNDO26WrVUSWA0Nc7wcDK1rJllJhY7Yxd1RszEUec2cXiLeB4SKyk4h0IZ7IeqKVbYoaG7ONuaNiYy7ymDuNQIjINOANYISILBeRc1S1CbgYeA6YD/xdVee1pp3FxMZsY8bGbGMu5J6JqVCGYRiGkUKn8SAMwzCM3DCBMAzDMHwxgTAMwzB8MYEwDMMwfDGBMAzDMHwxgTAMwzB8MYEw2gQioiJyq+v9ZSJybZH6/j8ROaEYfWW5z4kiMl9EXgp5/jMi0rdI996U5XhfEbnQ9X5gZ9r7w8gPEwijrVAPHC8i/VrbEDciUp7D6ecAF6rqwWFOVtXxqvplXoblTl8gKRCqukJVIxdNo31jAmG0FZqAu4BLvAe8HoDzaTmxa9YrIvJ3EVkkIjeKyAQR+Y+IzBWRnV3dHJLYYWuRJHYVE5FyEblFRN4WkZiI/I+r35dE5EFgro89pyT6fz9RLhwR+QXwX8CfROQWz/nbi8irIjI7cc2BifaPRaSfiAwVkQUiMiVx/AEROUREXhORxSKyT+L8a0XkMle/74vIUM+9eorICyLybsJGp5T7jcDOCRtuSdzz/cQ1VSJyb+L890Tk4ET7mSLymIg8m7DjZtfP7f8S958rImm/M6Nj0NmruRptizuAmPMgCskoYFfidfKXAlNUdR8R+RHxct4TE+cNBQ4CdgZeEpFdgDOADaq6t4h0BV4TkecT5+8D7K6qH7lvJiIDgZuAvYjX2n9eRI5T1etF5FvAZao6y2PjqcBzqvrrhEfS3WccuwAnAucTr69zKnHBOQb4KXBcyJ9HHfAdVd2Y8MbeFJEngCsT4xmdGMdQ1zUXAajqHiIyMjGm6sSx0cAY4h7eQhH5A7AtMEhVd0/01TekbUY7wzwIo82gqhuB+4Af5nDZ26q6UlXrie+s5Tzg5xIXBYe/q2qLqi4mLiQjgcOAM0RkNvAWsA0wPHH+f7zikGBv4GVVXZuog/MA8Z2+MtoInJXIqeyhqjU+53ykqnMTexXMA17QeB0c7ziyIcBvRCQG/Iv4fgEDslzzX8BfAVR1AfAJ4AjEC6q6QVXrgA+AHYn//IaJyB9E5AhgYw72Ge0IEwijrfF74rH8Hq62JhJ/qyIixLdWdKh3vW5xvW8h1UP2Fh1T4g/TH6jq6MTXTqrqCMzmAPv86u9nJLET2DeBz4C/isgZPqeFGUfy55CgyqefCUB/YK+Et7A64Dw3mcbktqsZqFDVL4h7bi8T9z6mZOnfaKeYQBhtClX9nPjm8ue4mj8mHtKB+BaSlXl0faKIlCXyEsOAhcQrYP6viFQCiEi1iPTI1AlxT+OgRO6gHDgFeCXTBSKyI7BGVe8GphLfVzgfPnauFZE9gZ18zumTuFdjIpewY6K9BugV0O+rJLYfTYSWdiD+8/ElEboqU9VHgZ+T/3iMNo7lIIy2yK3ESxg73A1MF5H/AC8Q/Ok+EwuJP8gHABeoap2ITCEevnk34ZmsJUusX1VXishVwEvEP3k/o6rTs9x7HPATEWkENhHPfeTDo3wVEnsbWORzzgPAkyIyC5gNLEjYvT6R9H4f+CfxfI/DH4kn1+cS91LOVNX6+I/El0HAvfLV3s5X5Tkeo41j5b4NwzAMXyzEZBiGYfhiAmEYhmH4YgJhGIZh+GICYRiGYfhiAmEYhmH4YgJhGIZh+GICYRiGYfhiAmEYhmH48v/VKeXgcPLhagAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(n, mean_list, color='black') #https://stackoverflow.com/questions/41709257/how-to-change-the-plot-line-color-from-blue-to-black#:~:text=The%20usual%20way%20to%20set,explicitely%20stating%20the%20color%20argument.\n",
    "#plt.xticks(range(1,iterations+1)) #https://stackoverflow.com/questions/52229875/how-to-force-matplotlib-to-show-values-on-x-axis-as-integers\n",
    "plt.xticks(n) #https://towardsdatascience.com/change-tick-frequency-matplotlib-axis-6af2c6bce1ea\n",
    "plt.xlabel(\"Number of simulations\")\n",
    "plt.ylabel(\"Unit weight sample mean\")\n",
    "\n",
    "plt.xscale('log') #https://stackoverflow.com/questions/773814/plot-logarithmic-axes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAEKCAYAAAAYd05sAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAA9GElEQVR4nO2deZgU5bX/P4dZWIdlkMVhCagscRfRYAQEgYCKokByXeISF1xRMfnlYnIlRp94YxK3XBJxi3G9xoWAM66giIDXBdS4hEUBUQGZYV+G2c/vj+5qqnt6qe7p6pnpOZ/n6YeuqrfeOu/M8K1Tp857XlFVDMMwjOyjVWMbYBiGYfiDCbxhGEaWYgJvGIaRpZjAG4ZhZCkm8IZhGFmKCbxhGEaWktvYBrg56KCDtF+/fo1thmEYRrNhxYoVW1W1W7RjTUrg+/Xrx/LlyxvbDMMwjGaDiGyIdcxCNIZhGFmKCbxhGEaWYgJvGIaRpZjAG4ZhZCkm8IZhGFmKCbxhGEaWYgLfwvnuu+/YsmVLY5thGIYPNKk8eCPzXHrppeTl5TF//vzGNsUwjDRjAt/CKSsro23bto1thmEYPmAhmhZORUUFNTU1jW2GYRg+YALfwqmoqKC2traxzTAMwwdM4Fs45sEbRvZiAt/CqaysNIE3jCzFBL6FYx68YWQvJvAtHBN4w8heTOBbMLW1tVRXV5vAG0aWYgLfgqmsrAQwgTeMLMUEvgVjAm8Y2Y2vAi8iM0TkcxH5TET+V0Ta+Hk9IzkqKioAE3jDyFZ8E3gR6QVcDwxV1SOBHOBcv66XrUybNo0333zTl74dgbeJToaRnfgdoskF2opILtAO2OTz9bKKuro6HnroIV5//XVf+jcP3jCyG98EXlU3An8CvgY2A7tUtZ5Sicg0EVkuIsvLysr8MqdZUl1dDUBVVZUv/ZvAG0Z242eIpgswCegPFAHtReSnke1U9UFVHaqqQ7t16+aXOc0S5yWo82+6MYE3jOzGzxDNWGC9qpapajUwF/ihj9fLOhxh98uDtywaw8hu/BT4r4FhItJORAQYA6z08XpZhyPsfnvwtbW1qKov1zAMo/HwMwb/HvA88CHwafBaD/p1vWzEbw/eEXiwTBrDyEZ8XdFJVX8D/MbPa2QzmYrBQyBMk5trC3wZRjZhM1mbMI7nngkP3uLwhpF9mMA3YTLpwVuIxjCyDxP4JkymsmjAPHjDyEZM4JswmcqiARN4w8hGTOCbMJnMojGBN4zswwS+CZPpLBrDMLILE/gmjGXRGIbREEzgmzDmwRuG0RBM4JswlkVjGEZD8CTwIjJcRH4W/N5NRPr7a5YBlkVjGEbDSCjwIvIb4D+Bm4O78oAn/TTKCGC1aAzDaAhePPhzgLOAfQCqugko8NMoI4B58IZhNAQvAl+lgVqyCiAi7f01yXBwhL2urs4XD9sE3jCyGy8C/6yIPAB0FpErgIXAQ/6aZUC45+5HmKaiooJWrQJ/AibwhpF9JKwPq6p/EpFxwG5gEDBLVRf4bpkRJuqVlZW0bds2rf1XVlbSoUMHdu/ebQJvGFlIQoEPZswscURdRNqKSD9V/cpv41o6mfDg27dvbwJvGFmKlxDNc0Cda7s2uM/wGbfA+/GitaKigg4dOgAWojGMbMSLwOeqash9DH7P988kw8HttfvlwZvAG0b24kXgy0TkLGdDRCYBW/0zyXAwD94wjIbgZRHOq4CnRGQ2IMA3wEW+WmUAmYnBOwJvE50MI/vwkkWzFhgmIh0AUdU9/ptlQP0smnTjZNGAefCGkY14yaJpDUwB+gG5IgKAqt7mq2WGrx58TU0NtbW1tG/fPrTt8Pzzz3P00UczcOBAANasWUOnTp3o0aNHWm0wDMNfvMTg5wOTgBoC5Qqcj+EzlZWV5OTkhL6nE2cWazQP/vLLL2fOnDmh7cmTJzNr1qy0Xt8wDP/xEoPvraoTfLfEqEdVVRUdO3Zkx44dnj34nTt3MmzYMJ566imOP/74mO3iCXxFRUXY9Xbs2MGePRaZM4zmhhcP/h0ROcp3S4x6VFZWUlBQEPruha+//prVq1ezYsWKuO0cgY8M0agqVVVVYYJfWVlpMXrDaIZ48eCHA5eIyHqgkkAmjarq0b5aZlBVVRUSeK8e/P79+4GAJx8PxyPv3LkzcEDga2trUdWwrJrKykrLsjGMZogXgT/NdyuMqFRWVtKtW7fQ92h89NFHrFixgssvvxzwLvC7du0CoGvXrsABga+urg7bdq5tAm8YzY+EIRpV3aCqG4D9BEoGh0oHG/7iDtHE8uAfeeQRZsyYEdpuqMA713GHbKqrqz0J/MUXX8z06dMTtjMMIzN4WdHpLBH5AlgPLAa+Al7x2S6D8BBNLA++oqKCffv2ESjZn7rAOwIeKfCR2/F4/PHHmT17dsJ20SgtLWXZsmUpnWsYRnS8vGS9HRgGrFHV/sAYwP4nZoDKyko6duwIxPbgKysrUdXQDSDdHrzTr98hmgkTJjB8+HB7mWsYacSLwFer6jaglYi0UtVFwLH+mmWAtywaJxtm377A1ITy8nIgfTF4R/D9Fvi1a9cCsHHjRl+vYxgtCS8CvzNYpuBtAjVp7iMw6cnwmaqqqlCeeiwP3hF4R9gdD37Hjh1x+961axetWrWiQ4cO5OTkNJoHP2nSJK688kr69OkDwFdffRW3/eeff862bdt8scUwsg0vWTSTgApgBnAB0AmwMgU+U1NTQ11dHW3btiUnJ8ezB+81RLNz5046duyIiJCbmxvTY/db4F988UUARo0aBUBZWVnc9kceeSS9evXi22+/9cUew8gmvBQbc5cleMxHWwwXjrC2bt2a1q1bJ+3BewnRODnwXjz4TMXGvdxILIxjGN6IGaIRkaXBf/eIyG7XZ4+I7M6ciS0TR2jz8/PJz89P2oOvqKgIHYvGrl276NSpE0CYB99YMXiHurq6xI0Mw/BETA9eVYcH/y3InDmGg1cP3mkX6cFDQMTbtGkT9bxYAp/JGPyUKVNC3500TxN4w0gfcV+yikgrEfksU8YYB3ALfCoePMQP0/gp8I5YJ2Lu3Ln1zrEZs4aRPuIKvKrWAf8Skb6pdC4inUXkeRFZJSIrReSklKxsgbhDNKnE4CE5gU800SkZ4XXCPNHYu3dvXC891jFVjRty8sqaNWu46aabPN+EDKM54yVN8mDgcxF5Q0RedD4e+78PeFVVBwPHACtTNbSl0VgefGQMPpWXrLGEePv27RQUFPC73/0u5rmxBP5Xv/oVbdu29WxDLCZNmsQ999zDmjVrGtyXYTR1vKRJ/jaVjkWkIzASuARAVauA9C8smqWkmkVTXl5Ou3btKC8vjynwqupriKaioiI0A9fN119/DQRWjLrlllvq2QSxBf6JJ57wfP14OONyViYzjGzGS7GxxdE+Hvo+BCgDHhWRj0TkYRFpH9lIRKaJyHIRWZ4oB7ol4TWLxtnv9uCLioqA2B58eXk5tbW1cQW+IXnwsTx4x0anBn00amtrWbp0ab3xOlU1HT744APP9kTDBN5oCXgpNjZMRD4Qkb0iUiUitR7TJHOBIcD9qnocgWX+ZkY2UtUHVXWoqg6N/E/ckvHiwdfV1YX2u2PwiQTeKVPgxYNPJQafisA7HvzKlSsZMWIE119/fdhxp6SCw4knnpgw1z8aFns3WhJeYvCzgfOAL4C2wOXBfYn4FvhWVd8Lbj9PQPANDzjCGi8G797n9uALCwvJz89PSeBjxeAzJfBbt24F4JNPPgk73qpV/T/VVBYid65jHrzREvAi8Kjql0COqtaq6qPAKA/nfAd8IyKDgrvGAP9O1dCWhiOs8bJo3ELq9uDbtm1L586dPQt8KjNZZ86cyY033hi1f7ddqsqoUaMoLi5m9+7Ag1+7du1ijPpADD7S047meTdEpE3gjZaAF4EvF5F84GMR+YOIzABiB1HDmU6gQNknBCpQ3pGamS0PL1k0biF1e/DJCnwqIZo777yT++67L7TtFmC3XZWVlSxevJipU6eGrhvPg/d7opN58EZLwovAXxhsdx2BOHofYErcM4Ko6sfB+PrRqnq2qsYvcWiE8JIHn8iDj1VR0o8QjVuY3Xa5c+Kd60bz4CMF3osH/9prrzF37txQqWEI3BBnzqz3qqdeP26B/+Uvf8nixV7yBgyjeeFF4IcQWGR7t6r+VlVvCoZsDB/x4sHHisGn4sHHmujkVeDdx925+I7Ai0jouvG853gTnSK58MILmTJlCsccc0xoX1VVFXfeeWfM/qP188c//jFUzTJdLFu2zMouGI2OF4E/C1gjIk+IyBki4iV33mggXrJoIkM0qpq2EE2yaZJeBN6xJ9pMV+f8VEoVODc3LzgCf9999zF//vyox9zcd999/Pvfyb06WrBgAcOHD+fuu+9O6jzDSDde8uB/BhwGPAecD6wVkYf9Nqyl4yUP3hH4/Px8ysvLQ23atWuXUOCdxT7AWww+0UxW93H3zcj57vbgo/XliH660hida8Xivvvu4+yzzw7bN2/ePCBwk3nhhRdQVW688UaOP/549u7dy5QpU1i9enW9G0plZSV79+4NbX/zzTcASd8YDCPdeM2iqSaw0PYzwAoCi4AYPpKMB9+1a1f27dsXisMn8uD37NlDhw4dQqGSdMTg3cedPr744oswu+MJvNPOHYOvqalh7ty5CUU/Wshn4sSJUc+L3HfEEUeEvm/evBmAe++9l6lTp4Zmz1ZUVFBQUMDcuXMZPHgwRx99dFgfPXr0CC2t6LbHcu6NxsbLRKcJIvJ34EtgKvAwgfo0ho8kk0VTWFhIeXl5KDTiCHxlZWXUnPTy8vKwTJZ0lCqIFPgVK1YwcOBA7rrrLiBxiCZS4CEQG58yZUpCkY9WEnnp0qVcccUVQEBonfMj+3F72Y4wOwuKXHzxxVGvt27dutD3VatWxXxaMIE3GhsvHvwlwDxgoKperKovq6qtyeozjuDl5ubSunVrqqur6wmGI76FhYXs27evnsBD9Nms+/btSyjwtbW1qGqY8MYTrEiBd9ZWXbBgAZA4RBMZ+1fVUO2a0tLSuNdu3bp11P2PPPIIEJgkNXr06FC/sXAE/vHHH4/Zxt12xYoVfPHFF/WOvfTSSwnPd9OvXz+uuuqqsH379++PWZ7CMLziJQZ/rqrOU1X7a8sglZWVtG7dGhEhPz8fqD9z0x2iifTgu3TpAsQWeHeqYjSBh4DYukUmXlaIW7Srq6tDouuOTTsCH82Dd9ZYdffjiPHXX3/NW2+9FfPajgcfT7ydNEgvXrXXRb2HDh3KpZdeGtpesWIFAC+88ILnawFs2LCBBx54IGxfu3bt6N+/f2h7/fr15ObmkpubG3fd2pqaGntyMEJ4isEbmccReDjgocYS+MLCQqqrq0MzRRN58JEhGvdMVrf41tTUhAl8vBetkR68c1NyBD6RB++cv3DhwnrHEmWjtG7dmi1btnD88cfHbZcobXHr1q1s2LAhbpto5zgMHTqUfv36hbaTFVpV5c033wydt3nz5tDP//e//z21tbXU1taGnoqikZeXV+9pAAI3rUytq2s0HUzgmyhVVVUhkXT+jXxkd3vwcMDzdAt8tMlOqXrw0eLw3bp1Y9WqVVFfssIBga+pqQnZ60Votm3bFhK6eKUNAHbv3k3Pnj356KOP4rbbs2dPXNG95ZZbwgQ6Fdw3iCeeeIJNmzZ5PnfOnDmMGTMmrO5O37592bhxI+vXrw/tc49h7ty5PPzww2E3mgcffDCs3+rqag466CCmTZuW1FiM5o8ngReRtq6aMkYGSNaDh3CBd7I63CESh2gvWSMnOkFAiCMFf+vWrcyZMye0b+vWrdx///31QjSxbkbO8USsX78+lBUUa11Zh1gzdqO1czJlMsU111wDBH6f48ePD92EysvL2bp1ayg0BfDcc8/VO7+0tJTevXuHee3uJ5EpU6ZwxRVXMHLkyJg2OL/Df/zjHw0bjNHs8JJFcybwMfBqcPvYJFZ0MlLELfCJPHhH4B0vrl27diGv1xFJN9E8+P379zN79uywHO/IEE1FRQXdunXj6quvDuuvrq6ungcfr9Kj11CBO+0zVdzrvl555ZUp95Mq8+fPp3///rRp04bXX3+dq666ipUrV9K+fXu6devG/fffH2q7aNEiT306Ar9s2bLQvpUrV8Z8d+C0Ly8v5/e//32qQzGaIV48+FuBE4GdEKgvA/TzyyAjgDtEE8uDd2fRwAGBb9u2bchDjzbLM5oHX1ZWxvTp03nzzTdD+yMF/tNPP41qa6TAV1VVxc0ASTYWnMiDj8eUKQfKJr3++usp99MQnIwiCPxsDj/88ND2HXckX39v/vz57N27l+HDh4ft/81vfhO2/eqrrzJ79uyw383NN9/MunXrwp7CHO68886Yv2OHl19+maKiorSsj2v4jxeBr1HV+NMCjbTj1YPPy8sLhWPcAp+sB+/gfvyPDNG4SxC4SdaDr66u9hQqcW4E2ZQV8vHHHze4j3/+859hE6sc3KUXjjjiCE477TSmT58eyqhyOPTQQ7n66qvZv38/l112GQMGDKC2tpaZM2dy9NFH89133/Hzn/+c8vJyVJWNGzfyzjvv8NVXX3HjjTeyefPmpF5G9+zZk3vuuSf1ARsp46WuzGcicj6QIyIDgOuBd/w1y6iqqvIUg2/Tpk1IrN0xeGefVw8+GpEefCyB37JlS9jsTi8hmsgyAbHaQWAykZEYdzzfS5mEd999l7/97W9AeGrowQcH5jHGy166+eabuemmm+jZsycDBgzgjTfeYPTo0YgIL774IjU1NUyePBkI/H3cdNNNbNq0icGDB/Ozn/0s6gIuRvrx8lOeDhwBVAL/C+wGbvTRJoOAt+4li6ZNmzYhsXZ78K1ataJNmzb1PHjHK4/lwUe29SLw7gwPiP6S1c3777/P+++/H/O4+/qGf5x66qmh7z169Ejq3H/+85+MGDGCAQMGADBmzBi6devGvHnzmDRpElOmTGHXrl1hs3z/9Kc/cfnll9O3b19EhFmzZjF//nx2797N119/HZo8ZqQPLxOdylX116p6QrC2+69V1QJwPuM1i6Z169ZRPXgIvGyN9OAdwU/Gg8/JyQFiC3xkLZhEHrxXTOCbF9u2beOcc84JbXfu3DmUruvGKQVx++23c/bZZ9OpUye+973vATBhwgRuuOEGysrKuPbaa8nNzUVEQplXX375JUuWLAn7uy4rK4uaLWbEEXgRKRaRF2N9MmlkS8RrHnykB5+fnx96/G3fvn09D975j+H24CNvAo5g19bWhnn76RD4ZFZS8pJOaWQXW7du5c9//jPdu3fnr3/9a+jdTn5+PiLCgAEDGDlyJB06dGDUqFGICN27d6egoID3338fEWHYsGFUV1ezdu1aVJW9e/fy/PPPh64xcuRIHn300cYaYkaJF4P/U8asMOrhxYOvrKwMi8Fv37497OWbVw/evSISBJ4AysvLQx58u3bt2LNnT0yBjxTiffv28fTTT0dtG82mWJjAG/GIXIXrBz/4AQDvvfdeyCmKxZIlS8LKTDj8/Oc/54EHHmDv3r1cdtllrFu3jn/84x+8+eabTJw4kZ07d7J06VLOOeechNdoCsQUeFUN/fSCa7IOBhRYraoNf/424uI1i8btwdfV1YXljHv14B2BP+SQQ1i3bh3t2rULE/iDDjoIiO3BR9oVuZCGm7Zt23oWeAvRGJnGqX4KB4rVde/ePaW+hg0bxhVXXEG/fv047rjj6mUzZQIvE53OANYCfwZmA1+KyGl+G9bS8ZIH7wi8U5QMwoXbqwd/7bXXAoE/SHcfTqjFuWl4Ffh4JCo74MYt8L179/Z8nmE0Bd59910uu+wyxowZQ2FhISIS8+PXJDwvWTR3AaNVdZSqngKMBiyp1WeS8eBFJCTYqXjw06dPR1Xp1atXWB+OoCeKwScj8M4qUl5w54z/+Mc/rne8oXVjDKOp8OCDD/oSkvQi8KURi2yvA0rTbkkLoLq6Omr98Ggkk0UDB0TYLfBuD/7+++9nxowZoW23B+/ghGIic+gTCbx7VmOiWafRJuh4IVrN91QfnQ2jqXHwwQeTl5eX9n69CPznIvKyiFwiIhcDxcAHIjJZRCan3aIs5umnn+bwww8PrdkZDy9ZNM5LViChB//KK6/wzDPPhLajhUqcqpROH5FtvXjwiQQ+GQ/eTbQ/fid90zCaMuPHj+eMM87giiuuYP369aHFdNyfZKqOJoOXmaxtgC3AKcHtMqAQOJPAS9e5Mc4zIli/fj01NTW8/fbbXHDBBXHbevXgHUFN5MHv3r07LF84ngfvp8Cn6sFHE/hkUi4NwwuHHXYYX375JdOnT2fnzp088cQT3HTTTUyaNIndu3dzwgkn0LFjR9q2bcv+/fupra2lQ4cO7N69m4qKiib3VJlQ4FX1Z5kwpCXgrMSzdOnSuAJfV1cXtipSohg8RPfgnWwYCAh8bW1t6OkhmgefKEQTq8CUu35NosqPqXrw0bx1Px5pjexg4MCBrFmzhkmTJvHZZ5+xePFiVq1aRWlpKeeddx6qyj/+8Q+mTp0ac6IfxF++0f233rFjRzp27JjWMaSDhAIvIv0JlCvo526vqmf5Z1Z2UloaeHWxZMmSuO2cly2OsDuz+ZL14Nu3b8++fftQVfbs2QMcKCsQzYN3QjTJxuDd+OXBR1uZKlZGzo9+9KNGqxxppM6ll14aqo3z5JNPMmzYML744gs6dOjAypUrGT16NIcddhgAzzzzDK+//jp5eXlMmzaNAQMG0LFjR6qrqxGRqKLtJBFA4Onv3HPPzczAGhEvIZp5wCMEYu/x1zwz4uJ48J9//jnbtm0LCWokjqfuePDOuqypePB1dXVUVVWFlvNbv349IhJViHv06EGrVq1C5YcjQzTl5eVhqz9FoyEx+LPPPpt58+ZFPRYtdz6WwE+fPt0EvhGYOHEiJSUloe0BAwbQo0ePuA7Nrl27eOqpp7j66qsRESZPnsy3334besI99NBDAeqVRj733HOjCrQ91YXj5SVrhar+WVUXqepi5+O7ZVlIWVlZqKjTO+/ELsgZKfDOdy9ZNG7Rc0S/vLw8TODbtWsXNX7dpUsXli5dyiWXXAJE9+AT/QdyC7xzs3HP+HME3l0T3SHeqkTRnh4ashCI0TB++9vfhm2rKsXFxdxxxx0cddRR/Otf/2LNmjUJn1Y7derENddcE/p7POOMMxplYZZsxYvA3ycivxGRk0RkiPPx3bIspKysjAkTJpCfnx/3D98RcrcwRnrwNTU11NbWJvTgIeAlOfHzTZs2xZ1sdNJJJ9GpUycgegw+0fRs9/WdxUPcNwXH3mjlYp13ANGI5sGbwEfn/PPPD33v3bs3U6dODTv+2GOPUV5ezpo1a+rV2j/11FP57rvvQtuVlZVRlxK85ZZbePrpp9m+fTtbtmwJ7b/55pv55JNPwspHG42HF4E/CrgC+D2BSU93YXVqksZZz7RPnz6ccMIJcQXeiwfvtEkUgwfC/sOqatT4uxsnfhm5ZN7+/fsTCrxbwB2b3Oc4Y4r20jRWyMptixv3ak0Ao0aN4qGHHor6dHLcccfFtTubePzxx9m8eTPV1dV88803PPfcc2F1Wy688ELatm0bKvXr5uSTT6ZHjx6h9L38/HymTp1KcXFxqM2VV16JiHDeeefRpUuXJpc5YhzAi8CfAxyiqqeo6ujg59SEZxlhbN++HVWle/fuDB8+nBUrVkQVLTjgwbsFPtKDdzxyLx585OpJicoFRAp8MiEa5/q5ubkh29znOGIfzYN3Yv9w4H2Fg7t2OQTeY4wfPz5s35NPPsnll18e1a5nn32Wbt26xbW9oTz88MP19r366qu+XtPBHTLJycmhZ8+eYS8a3eGvWOmlq1atqrfsn8PEiRPZtm0bt912G3/961/TZLXhN14E/l9AZ5/tyHocwerWrRsjRoyguro65qIXjpBHer5uDz5S4L168O79sXC862gx+EQevHNubm5uqK1b4N0vjiNx39Aiwy833nhj1Ou4iZcX3717dxYuXBjX9oZy2WWXsXXrVh577DEWL17MwoULYxaYeuyxxwCYNGkSAKeffnrcvlWV5cuXx7w5O08oJ510Usw+Tj755LjXGDRoUNzJY4WFhdxyyy22GlMzwstvqgewSkRes3rwqeMW+B/+8IeICEuXLo3aNlqIpiEevCPwjkfnNUQTWdbAi8A7Yp6XlxeK77ozZ5wxRRMJt7hEPilEinc0MY8nPK1atUqrMB1zzDFh2076Z9euXbnooosYOXIkY8aMiXquqnLRRRehqvztb3+jbdu2zJw5s167WbNmhW0ff/zx7Nu3j//6r/+q13bChAmsXr2a1157LabNr7zyCp9//nm9/bfcckvcXHCj+eLlL/43BMI0d3AgBn9X3DOMejg58N26daNLly4ceeSRMePw0V6yxorBJ6pFAwcEvn///mH7YxErBu/EZGPxxBNPhAqA5eTk0K9fP375y1/y4osvhkTJEW5HzLt06cK7777Lc889F1fgI4m2ELcj+rHEP9XSBsuWLau3b9GiRZ7OTbRgeGFhIeXl5YwYMYLVq1eH9t96663ceuutUc+5/fbbqa2tZcaMGaxYsYKamhry8vIYOHBg3HkGBQUFUbOXbrvtNqu9n6V4WbJvcbRPJozLJhwP3nkhNWLECN55552oOeXp8uCdfU4M3nmplqwH774hxBLe1q1b89Of/jRkj1MG9c4772TgwIHcfvvtqGpotp+zsHNdXR0/+MEPmDp1apgAu0X6kEMOiWtv5DnRBD4nJydlD/6HP/xh2PY333xTL/QSS8iPOuooDjnkEM4444yE1xk4cCAbNmzgtttuY9asWYgIp5xyStS2rVq14u6772bIkCFWk8eIiZd68MNE5AMR2SsiVSJSKyK7M2FcNuEIvJMpMmLECPbu3csnn3xSr20sgfcSg4+sBw8HPHhnFqBXDz6awMfy4B3hd+xxllqLZPz48cyZM4d77glUnHaXOYglVLFCWZE4Aj5o0CAA5syZE3YsnhAmU3o4Wm36aCEWCPzs1q5dW+8lcSz69u3LLbfcErpJvfzyy2zYsMGzbYbhxotLMxs4D/gCaAtcHtxnJEFpaSldunQJCaEzMy9amCZWiCZVDz5S4BviwScSeOemFGu2q7O4gZNr70XgHW/fjXvaubtvCHj8e/bsYdq0aaFjiWLwI0aMCNu+++67Y7aNRFX59a9/HbdNqsu7tWvXjr59+6Z0rmF4emYN1oPPUdVaVX0UGOX1AiKSIyIfiUhJ4tbZS1lZWVi+cO/evenXr19UgU/Fgz/66KMZMmRIWIw10oN3pn2nmiYJASGfN29eaEaue7/bnlgefOQ13ALvtU6NqkYteeAW8A4dOiAinHXWWaFjsW4gTokGNzNmzACoJ65///vfPdkYSbR69obhN14Evjy4JuvHIvIHEZkBxHcBw7kBWJmSdVlEWVlZvTzsESNGsHTp0nrxWy8TnSIFvk+fPqxYsYKePXuG2uTn55OTk0N1dTUFBQUhT9hrmmTkS1anz0mTJvHii+GJVI6H2hCBj5annkxtkWix9+eee46ysjJEJKoHf8opp7Bs2TJuvPHGeje+VatW8dlnn4Xtu/jiiz3b4+bCCy/k/PPPr7fAuWH4iReBvzDY7jpgH9AHmBL3jCAi0hs4A6g/A6SFEU3ghw8fzpYtW/jyyy/D9nspVRDtJhCJiIREq6CgIPQEkciDFxFycnJQVXJzc6NOVIr0hh1xTRSicXDOj7y5XXzxxaGMmxUrVoSqX3ohmoDn5+eHSiBEO37ddddx6KGHcuyxx9YrhzBo0KDQU8WiRYs8r8YVjTZt2vDUU095fmFsGOnASz145w1PhYj8GegTsYRfPO4FfgnEfPYWkWnANKj/OJxNlJaW1pto4sR9lyxZEjZtPBUPPhbt27dnz549dOzYkaKiImbOnBkKW8QjJyeH2tpaWrduHZYjHUvgnXK+jj1uzzwaeXl5XHnllVx00UVh+90hkCFDkit5lGgBkGghGq/1bEaNGpWULYbRFPCSRfOWiHQUkUICs1ofFZGEb6BEZCKB9VxXxGunqg+q6lBVHer3VPLGoq6ujm3bttWr2TF48GC6du1aL0sklTTJWDjeeseOHRER/vu//5vBgwcntNkRdSfM4+B485ETY5x684nscRAR5syZUy8FsSEkEvhoHvxpp50Wtn3//ffHrfRpGM0JL9PXOqnqbhG5HHhUVX8jIvVz++pzMnCWiJxOYNm/jiLypKr+tCEGN0e2b99OXV1dvRCNiDB8+PB6Au9lolMyHjyQ9GozjoC3bt06TOBjefBuOxuLRHnukTYfddRR9c656qqr0m6XYTQWXmLwuSJyMPATwHMmjKrerKq9VbUfcC7wZksUdwifxRrJ4YcfzldffRUW0vDDg092JaVUBd6rB+8HyXrwtqarke14EfjbgNeAL1X1AxE5hEBOvOGRyFmsboqKiqiurmbr1q2hfZWVleTm5oYJUqxSBYnyq90hmmSIFaJxLyMYjaYs8LFeDBtGtuKlVMFzqnq0ql4T3F6nqp6yaFx9vKWqE1M1srnjLjQWSVFRERBYiMOhqqqqnnDn5+dTV1cXyk5xVnNKJFLpDtFE1pKJJNMhGmf1KS+YB2+0NKzuZwZIVuArKyvrCaWz7Xjx7vVY49FQD965iTji2Jghmk2bNtVLm3zkkUdCtdAT5cybB2+0NEzgM4ATg4+2JF0sDz5S4B1hdUIzXgXe8eCTjcFHlh5wxDGWwDuzZP0U+IMPPrhezZhWrVoxa9YsVDVh0a1ID97qmhvZTsIsGhHpr6rrE+0zYlNWVhZWh8aNM7s00oOPDNFk2oM//vjjWbt2bT2bo6VJ7t+/PySWjRmDT0S6PfiPP/446lqxhtFU8OLCvBBl3/PpNiSbiTaL1SEvL4/u3bsnDNE01INPVuDPPPNMAD799FOAUL3waB58mzZtQvubcs2VdHvsxxxzTFrz+A0j3cT04EVkMHAE0ElEJrsOdSSQ1254JJ7AQyBMk+gla6QHX1lZ6asH70wA2r59OwDHHnss3333Heeddx7QdF6yJkOkwCdaJs8wmjvxQjSDgIkE1mM907V/D3CFjzZlHaWlpQwcODDm8aKiIjZu3Bja9urBexHTVGPwXbt25fzzz+fEE08E4KOPPgo7HitNsrnEtTdu3BhWmM0wspGYAq+q84H5InKSqv5fBm3KOsrKyuIueFxUVMSHH34Y2m4KWTQATz31VMxjzX0VIefltmFkM15KFXwpIr8C+rnbq+qlfhmVTdTV1bF169aEIZotW7ZQU1NDbm5uzDx4CPfgvYh2YWEhED2DpyEkEngvS9Q1FrGWwTOMbMOLwM8HlgALgfhFvo16xKpD46aoqAhVZcuWLfTq1YvKysrQikcO0Tz4aDNjI5k8eTJvvPFG2it1xhP4/fv3J1XHPZNs2LAh7Tc7w2iqeBH4dqr6n75bkqXEK1Pg4M6FdwQ+XVk0+fn5ntcDTYZ4sfamnCqZzSWpDSMSL2/ESoIVIY0UiDeL1SFyslM6s2gMw2i5xEuT3AMoIMCvRKQSqA5uq6om/9auBeJF4J0FpB2BT2cWjd+cc845jW2CYRgxiJdFk1xenREVp0xBvBBNt27dyMnJiSvwqWbR+MnOnTsTLv9nGEbj4aVUQbR103YBG1Q1/sKbRsiDj/diLycnh549e4Zy4b1m0TS2wEe+CDYMo2nh5SXrX4EhwKfB7aMILN3XVUSuUtXX/TIuGygrK6Nz584Js0rcs1kTefCq2iQE3jCMpo2Xl6xfAcep6vGqejxwLPAZMBb4g3+mZQelpaVx4+8OiQTe7cHX1NSgqibwhmHExYvAD1bVz50NVf03AcFf559Z2UNZWZmnfHW3wCfKovG6XJ9hGC0bLyGa1SJyP/BMcPs/gDUi0ppAVo0Rh7KyMg477LCE7YqKiti2bRv79u2jrq4urgfv3AjsBadhGPHw4sFfAnwJ3AjMANYF91UDo32yK2tIVEnSwcmF37BhA1C/KqMj8FVVVfzlL38hLy+PiRNb7CqIhmF4IKEHr6r7gbuCn0j2pt2iLMKpQ+M1RAOElqSLDNG0atWK3NxcNm7cyJNPPskFF1xA796902+0YRhZQ7yJTs+q6k9E5FMCE57CUNWjfbUsC9ixYwe1tbWePHhnspMj8NEmMbVu3Zonn3yS/fv384tf/CK9xhqGkXXE8+BvCP5rcYAU8TKL1SHSg48m8Pn5+ezYsYOJEydyxBFHpNFSwzCykZgxeFXdHPx3Q3DXgOD3UmB7Bmxr9iQj8IWFheTn58cM0cAB0f/lL3+ZRisNw8hWvMxkvQKYBhQChwK9gTnAGH9Nyzx1dXWoatoWs/BSpsBBRCgqKorrwRcUFNCvXz+GDx+eFvsMw8huvGTRXAucDOwGUNUvgMSK1Qy56aabGDZsWNr6S8aDh0CY5quvvgKie/DPPPMMzz77LCKSNhsNw8hevAh8papWORsikkuUl67Nnbq6Op5++mmWL18eSlVsKF7q0LgpKipi586dQHQPfsiQIfTp0yctthmGkf14EfjFwSX72orIOOA5oNhfszLP+++/HxLkhQsXpqXP0tJSOnXqFNUbj4Z7ndCmUArYMIzmjReBnwmUESg2diXwMvBffhrVGJSUlJCTk0PXrl3TJvBeyxQ4uAXe603BMAwjFl5KFYwCnlLVh3y2pVEpKSnh5JNPpm/fvrz66qvU1dXFXZbOC15nsTqYB28YRjrxWqrgYxH5PxH5g4icKSJdfLYro3z99df861//YuLEiYwbN46tW7fyySefNLjfZAXemewEJvCGYTSchAKvqhep6kBgCvAt8BcCIZus4aWXXgLgzDPPZMyYQPbnggULGtxvaWmphWgMw2g0Egq8iPxURB4AnidQA342MMJvwzJJSUkJhx56KIMGDaJXr14cfvjhDY7DO3VoLERjGEZj4SVEcy+BRT4eAq5X1T+o6v/5aVQm2bdvH2+88QYTJ04M5ZePHTuWJUuWhOqup8LOnTs916FxKCgooH379oAJvGEYDcdLiOYg4FKgDfA7EXlfRJ7w3bIM8cYbb1BZWcmZZ54Z2jdu3Dj279/PO++8k3K/zizWZATemc0KFqIxDKPheAnRdAT6At8D+gGdgDp/zcocJSUlFBQUMGLEgajTKaecQk5OToPCNE5OfTIxeDgQpjEP3jCMhuIlRLMUOBP4BPgPVR2kqhf7a1ZmUFVKSkoYP358mMdcUFDAsGHDGvSiNdkyBQ4m8IZhpAsvIZqjVfUaVX1aVb/NhFGZ4sMPP2Tz5s1h4RmHcePGsWLFCrZvT61wZqoC37t3b3JycsjN9TJFwTAMIzYNm8kTBxHpIyKLRGSliHwuIjckPiuzlJSUICKcdtpp9Y6NHTsWVWXRokUp9Z1KDB5g+vTpVlDMMIy04JvAAzXAz1X1+8Aw4FoROdzH6yVNSUkJw4YNiyrCJ554IgUFBSmHacrKypKqQ+PQp08fJk+enNI1DcMw3Hh5yXqyl32RqOpmVf0w+H0PsBLoFf+szLFp0yaWL18eNTwDkJeXx6hRo1J+0ZrsLFbDMIx048WD/x+P+2IiIv2A44D3ohybJiLLRWS5E7fOBC+//DIAEyfGXpFw3LhxrF27NrQIRzKYwBuG0djEFHgROUlEfg50E5GbXJ9bAc9LHolIB+AF4EZV3R15XFUfVNWhqjo0k4JYUlJC3759OfLII2O2GTt2LJBa+eBkyxQYhmGkm3gefD7QgUDFyQLXZzcw1UvnIpJHQNyfUtW5DTM1fVRUVLBgwQLOPPPMuC8zBw8eTK9evVISePPgDcNobGLm4qnqYgKLffzdtfC2ZySgnI8AK1X17gbYmHYWLVpEeXl53PAMBGaWjh07lpKSkqTKB6dSh8YwDCPdxAvR3Bv8OltEXoz8eOj7ZOBC4FQR+Tj4OT0NNjeYkpIS2rdvz6hRoxK2HTt2LNu2bePjjz/23P/OnTupqamxEI1hGI1KvNk0Tr2ZP6XSsaouBZpcMreqUlxczLhx42jTpk3C9k4cfsGCBQwZMsTTNVKd5GQYhpFOYnrwqroi+O/iaJ/MmZhePv30U7755puE4RmHnj17cuSRRyYVhzeBNwyjKeApD15EFojIGhFZJyLrRWRdJozzg5KSEgBOP917tGjcuHEsWbKE/fv3e2rvhHP69OmTtH2GYRjpwstbw0eAu4HhwAnA0OC/zZLi4mJOOOEEDj74YM/njB07lsrKSpYtW5awbV1dHbNnz+aEE05g8ODBDTHVMAyjQXgR+F2q+oqqlqrqNufju2U+UFpaynvvvec5POMwcuRI8vLyPIVpXn31VVavXs2MGTOsnoxhGI2KF4FfJCJ/DE58GuJ8fLfMB1555RVUNWmB79ChAyeddJKnujT33nsvRUVFTJ3qaaqAYRiGb3ipSfuD4L9DXfsUODX95vhLcXExRUVFHHfccUmfO27cOGbNmsXWrVs56KCDorb57LPPWLBgAXfccQd5eXkNNdcwDKNBeKkHPzrKp9mJe1VVFa+99lrY2qvJ4JQPfvPNN2O2ue+++2jbti3Tpk1riKmGYRhpIaYHLyI3RexSYCuwVFWTr77VyLz99tvs3bs36fCMw9ChQ+nUqRMLFy7kJz/5Sb3jZWVlPPHEE1xyySV07dq1oeYahmE0mHgefEHEpyOBMM0rInJuBmxLK8XFxbRp04YxY8akdH5ubi6jR4+O+aL1gQceoLKykuuvv74hZhqGYaSNeLVofhttv4gUAguBZ/wyKt04s1fHjBlDu3btUu5n7NixzJs3j7Vr13LooYeG9ldVVfGXv/yF8ePHc/jhTWpNE8MwWjBJr+ikqttpgiUI4rFq1SrWr1+fcnjGYdy4cUD98sHPPvss3333HTNmzGhQ/4ZhGOkkaYEXkVOBHT7Y4hvFxcVA/MU9vDBgwAD69OkTli6pqtxzzz18//vf50c/+lGD+jcMw0gn8V6yfkrgxaqbQmATcJGfRqWbkpISjj32WHr37t2gfkSEcePG8c9//pPa2lpycnJYunQpH374IXPmzLGJTYZhNCniefATgTNdn4nAIFU9UVVXZcK4dLB9+3aWLVvWYO/dYezYsezYsYOPPvoIgHvuuYfCwkIuvPDCtPRvGIaRLuK9ZE16kY+myCuvvEJdXV3MxbWTxcnCWbBgAYWFhcybN4+ZM2c26OWtYRiGH3iZydqsKSkpoXv37gwdOjRxYw90796dY445hoULF1JWVkZOTg7XXnttWvo2DMNIJ0m/ZG1OVFdX8+qrr3LGGWd4Xm7PC2PHjmXp0qU8/PDD/OQnP6FXr15p69swDCNdZLXAL1u2jJ07d6YtPOMwbtw4qqqq2LNnDzfccENa+zYMw0gXWR2iKSkpIT8/P7TsXroYMWIE+fn5DB06lBNPPDGtfRuGYaSLrBf4UaNGUVBQkNZ+27Vrx/PPP8+AAQPS2q9hGEY6yVqB/+KLL1i9ejXXXXedL/2nO+xjGIaRbrI2Bu+svXrGGWc0siWGYRiNQ1YL/BFHHEH//v0b2xTDMIxGISsFfteuXbz99tsWRjEMo0WTlQL/2muvUVNTk7byBIZhGM2RrBT4kpISunbtyrBhwxrbFMMwjEYj6wS+traWl19+mdNPP52cnJzGNscwDKPRyDqBf/fdd9m2bZuFZwzDaPFkncCXlJSQm5vL+PHjG9sUwzCMRiXrBL64uJiRI0fSqVOnxjbFMAyjUckqgV+/fj2ff/65hWcMwzDIMoF/6aWXgIavvWoYhpENZJXAFxcXM2jQICsCZhiGQRYJ/J49e3jrrbfMezcMwwiSNQK/cOFCqqqqTOANwzCCZI3AFxcX07lzZ04++eTGNsUwDKNJkBUCX1dXx0svvcSECRPIy8trbHMMwzCaBL4KvIhMEJHVIvKliMz06zrLly+ntLTUwjOGYRgufBN4EckB/gKcBhwOnCcih/txreLiYlq1asVpp53mR/eGYRjNEj89+BOBL1V1napWAc8Ak/y4UElJCSeffDKFhYV+dG8YhtEs8VPgewHfuLa/De4LQ0SmichyEVleVlaW9EXKy8sBWyPVMAwjEj8X3ZYo+7TeDtUHgQcBhg4dWu94Itq1a8dHH31EXV1d8hYahmFkMX568N8CfVzbvYFNfl2sVausSAgyDMNIG36q4gfAABHpLyL5wLnAiz5ezzAMw3DhW4hGVWtE5DrgNSAH+Juqfu7X9QzDMIxw/IzBo6ovAy/7eQ3DMAwjOha4NgzDyFJM4A3DMLIUE3jDMIwsxQTeMAwjSxHVpOcW+YaIlAE7gV2u3Z3ibLu/HwRsTYMZkddLtW2sY9H2xxtj5LaNuWWNOV3jjWVTKu3SNWa/f8exbEqlXVMe8/dUtVvUI6rapD7Ag163I74v9+P6qbaNdSzafhuzjTnWmNM13mTGnKhdusbs9++4pY7Z/WmKIZriJLYjj/lx/VTbxjoWbb+N2cYcud2YY07ULl1j9nu8yfSbTWMO0aRCNA1BRJar6tDGtiOT2Jizn5Y2XrAxp5Om6MGnyoONbUAjYGPOflraeMHGnDayxoM3DMMwwskmD94wDMNwYQJvGIaRpZjAG4ZhZClZKfAi0l5EHhORh0Tkgsa2JxOIyCEi8oiIPN/YtmQKETk7+DueLyI/amx7MoGIfF9E5ojI8yJydWPbkymC/6dXiMjExrYlE4jIKBFZEvxdj0q1n2Yj8CLyNxEpFZHPIvZPEJHVIvKliMwM7p4MPK+qVwBnZdzYNJHMmDWwuPlljWNp+khyzPOCv+NLgP9oBHPTQpJjXqmqVwE/AZptKmGS/58B/hN4NrNWppckx6zAXqANgdXxUsOP2VN+fICRwBDgM9e+HGAtcAiQD/wLOBy4GTg22ObpxrY9E2N2HX++se1uhDHfBQxpbNszNWYCTss7wPmNbXsmxgyMJbAi3CXAxMa2PUNjbhU83gN4KtVrNhsPXlXfBrZH7D4R+FID3msV8AwwicAdr3ewTbMZYyRJjjkrSGbMEuBO4BVV/TDTtqaLZH/Pqvqiqv4QaLbhxyTHPBoYBpwPXCEizfL/dDJjVtW64PEdQOtUr+nrik4ZoBfwjWv7W+AHwJ+B2SJyBhmYDpxhoo5ZRLoCvwOOE5GbVfW/G8U6f4j1e55OwLvrJCKHqeqcxjDOJ2L9nkcRCEG2JvtWS4s6ZlW9DkBELgG2usQvG4j1e54MjAc6A7NT7by5C7xE2aequg/4WaaNyRCxxrwNuCrTxmSIWGP+M4GbeTYSa8xvAW9l1pSMEXXMoS+qf8+cKRkj1u95LjC3oZ03y0cdF98CfVzbvYFNjWRLprAx25izFRtzmsfc3AX+A2CAiPQXkXwCL2JebGSb/MbGbGPOVmzM6R5zY79ZTuIN9P8Cm4FqAne9y4L7TwfWEHgT/evGttPGbGO2MduYm8qYrdiYYRhGltLcQzSGYRhGDEzgDcMwshQTeMMwjCzFBN4wDCNLMYE3DMPIUkzgDcMwshQTeCMtiIiKyF2u7V+IyK1p6vvvIjI1HX0luM6PRWSliCzy2P5lEemcpmvvTXC8s4hc49ouakm1/43UMIE30kUlMFlEDmpsQ9yISE4SzS8DrlHV0V4aq+rpqrozJcOSpzMQEnhV3aSqvt/0jOaNCbyRLmqAB4EZkQciPXDHWw2uWrNYRJ4VkTUi8nsRuUBE3heRT0XkUFc3Y4Mr3KxxVvURkRwR+aOIfCAin4jIla5+F4nI08CnUew5L9j/Z8Fyw4jILGA4MEdE/hjR/mAReVtEPg6eMyK4/ysROUhE+onIKhF5OHj8KREZKyLLROQLETkx2P5WEfmFq9/PRKRfxLU6iMgbIvJh0EanFPTvgUODNvwxeM3Pgue0EZFHg+0/EpHRwf2XiMhcEXk1aMcfXD+3vwev/6mI1PudGdlBc68maTQt/gJ84giJR44Bvk+gTvY64GFVPVFEbiBQDvjGYLt+wCnAocAiETkMuAjYpaoniEhrYJmIvB5sfyJwpKqud19MRIqAO4HjCdTafl1EzlbV20TkVOAXqro8wsbzgddU9XfBJ4J2UcZxGPBjYBqB+iLnE7hhnAX8Cjjb48+jAjhHVXcHn4beFZEXgZnB8RwbHEc/1znXAqjqUSIyODimgcFjxwLHEXjCWi0i/wN0B3qp6pHBvjp7tM1oZpgHb6QNVd0NPA5cn8RpH6jqZlWtJFCLwxHoTwmIusOzqlqnql8QuBEMBn4EXCQiHwPvAV2BAcH270eKe5ATgLdUtUxVa4CnCKy0E9dG4GfBdwpHqeqeKG3Wq+qnGqhV/jnwhgbqgESOIxEC3CEinwALCdQL75HgnOHAEwCqugrYADgC/4aq7lLVCuDfwPcI/PwOEZH/EZEJwO4k7DOaESbwRrq5l0Asu71rXw3BvzUREQJLkzlUur7XubbrCH/CjCyapATEcLqqHhv89FdV5waxL4Z90epvx0UDK/GMBDYCT4jIRVGaeRlH6OcQpE2Ufi4AugHHB731LTHauYk3JrddtUCuqu4g8OT0FgHv/+EE/RvNFBN4I62o6nYCiyO7FwD/ikBIBAJLsOWl0PWPRaRVMC5/CLAaeA24WkTyAERkoIi0j9cJAU//lGDsPAc4D1gc7wQR+R5QqqoPAY8QWFczFb5yzhWRIUD/KG06Ba9VHYylfy+4fw9QEKPftwku3xcMzfQl8POJSjD000pVXwBuIfXxGE0ci8EbfnAXcJ1r+yFgvoi8D7xBbO86HqsJCHEP4CpVrRCRhwmEPz4MPhmUkSDWraqbReRmYBEBz/dlVZ2f4NqjgP8nItUEVrqP5sF74QUOhJQ+IFAiNpKngGIRWQ58DKwK2r0t+NL2M+AVAu87HP5K4OXwpwSeEi5R1crAjyQqvYBH5cDapjenOB6jiWPlgg3DMLIUC9EYhmFkKSbwhmEYWYoJvGEYRpZiAm8YhpGlmMAbhmFkKSbwhmEYWYoJvGEYRpZiAm8YhpGl/H/l/aouMD4XjgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(n, var_list, color='black')\n",
    "plt.xticks(range(1,simulations_number+1))\n",
    "plt.xlabel(\"Number of simulations\")\n",
    "plt.ylabel(\"Unit weight sample variance\")\n",
    "\n",
    "plt.xscale('log')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 102,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Probability density')"
      ]
     },
     "execution_count": 102,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEMCAYAAADTfFGvAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAA6PElEQVR4nO3dd3hc5Zn38e896lax5F5xw4BNIGAMgZACAQfjoDiQ3SykQNoSB1jSyIbdJBvS3iUJhCQbgkMNaRA6NphiG1NMc++yLVmWbXXJVq+jmfv94xyRYRhJM7Jmjsr9ua65Zua0+c14rHvOc855HlFVjDHGmGj5vA5gjDFmaLHCYYwxJiZWOIwxxsTECocxxpiYWOEwxhgTEyscxhhjYpLsdYBEGDdunM6cOdPrGMYYM6Rs3ry5VlXHh08fEYVj5syZbNq0yesYxhgzpIjIoUjTranKGGNMTKxwGGOMiYkVDmOMMTGxwmGMMSYmVjiMMcbExAqHMcaYmFjhMMYYExMrHMYYY2JihcMYD+Tn55Ofn+91DGP6xQqHMcaYmFjhMMYYExMrHMZ4zJqszFAzIjo5NGYwaOsM8OzOCt44UMuBmZ9EgO89toOG7Fl0dAVIS07yOqIxUbHCYUycBYLK394+xO0v7qehzc+4rFQEUITnd1fScMIlfOy2V/jepaeQf/pkRMTryMb0ypqqjImjupZOvnDf2/zP07s5bepoHvnaeWz8/sXMLlnBnJKn2fyDi5lxaBVjMlO58aGtfPMf22j3B7yObUyvbI/DmDgpr2/jqnveoqKhnall66jfvY+fPv3uZZKTfOQ0H+ap68/nD+uK+PWa/ZTXt3HvNWczOiPFm+DG9MH2OIyJg+rGdj5379sca+7koX8/lzH1++itASrJJ7z4m2/z+6sWsP1IA1/+00bb8zCDVlwLh4gsFpF9IlIkIjdHmH+KiLwpIh0iclPI9JNFZFvIrVFEvunOu0VEykLmLYnnezAmGqEX9LX7A/z7XzZT1djOn758NmfNyIt6O584fTK/ufIMthyu48aHthIMarwiG9NvcWuqEpEk4E5gEVAKbBSRFaq6J2SxY8CNwKdC11XVfcAZIdspA54MWeQOVb0tXtmN6a/L8vMpm/JR6vLmsfzzZ3HWjDExb2PJaZP5wSfm89Nn9vDHV4v5+gVz4pDUmP6L5x7HOUCRqharaifwMLA0dAFVrVbVjYC/l+1cBBxQ1Yhj3xozmDTkzKEubx43XHgii983qd/b+fL5M/nEaZO57cV9bCo5NoAJjTl+8SwcU4EjIc9L3WmxuhJ4KGzaDSKyQ0TuF5GI7QAicq2IbBKRTTU1Nf14WWNi40/OoHzyh8loreLl33/3uC7sExH+99OnMSU3nZse3W7HO8ygEs+zqiIdC4ypwVZEUoFPAv8VMvku4Kfutn4K3A58+T0vpHo3cDfAwoULraHYxF355I8Q9CUzrWwd4n7V+1M8QtdJzpzKwZn5/HZtId9bfMqAZTXmeMRzj6MUmB7yfBpQHuM2LgW2qGpV9wRVrVLVgKoGgXtwmsSM8VRj1gk05sxiYs0m0jvrB2y7WS1l5NUVcNe6QgoqGgdsu8Ycj3gWjo3AXBGZ5e45XAmsiHEbVxHWTCUik0OeXg7sOq6Uxhynzq4glZPOI7WjnrFHdwz49idVvUVSoJOfP1uAqu08G+/FralKVbtE5AbgBSAJuF9Vd4vIMnf+chGZBGwCcoCge8rtfFVtFJFROGdkfS1s078UkTNwmqpKIsw3JqH+/GYJHWl5zDi0Cp8GB3z7yYEOJtRsYn1yOi/vq+HCUyYM+GsYE4u4XjmuqquAVWHTloc8rsRpwoq0biswNsL0LwxwTGP6rbHdz+/WFpLVdJjs5sNxe50xdXtIP+0Sfr6qgNu/+1UEZeXKlXF7PWN6Y1eOG3McHlhfQmN7F5Oq3+71yvDj5dMg31t8CkXVzdSPnhvHVzKmb1Y4jOmnhjY/964v5uPzJ5LRfjTur/f7//p30ttqqR6/AI1rmTKmd1Y4jOmn+9YfpKm9i29efFJCXk+ACTWb6UzLpX70iQl5TWMiscJhTD80tft54PWDLD51EvOn5CTsdXOaDpLefpSa8QsIWD9WxiNWOIzph39sPEJTexeFK+5M6OsKML5mCx1peazeU9Xn8sbEgxUOY2LUFQjywOslnDNrDKPaE9+dzejGYlI6m7h//cGEv7YxYIXDmJg9t6uSsvo2vvqhWZ68vqCMPbaTDSXH2Fna4EkGM7JZ4TAmBqrKva8VM3PsKC6eN9GzHGPq9pKZmsR964vfNRaIMYlghcOYGGw5XMf20ga+8qFZ+HzenRKbFOzkM2dP55kdFfiTMz3LYUYmKxzGxOCvbx0mKy2ZKxZE7PAgob70wVkEVDmWN8/rKGaEscJhTJTqWjp5dmcFnzpzCplpce2tJyonjB3Fh+eO51jeKXZBoEkoKxzGRGnJ9T+msyvIZ8+Z4XWUd3z2nBPoSsmiKesEr6OYEcT7n03GDHL5+fkocOzEfyOjtSqhF/z1xsnlI/mkz3FszHyv45gRxPY4jIlCy6jJdKTlMbZuj9dR3kUIkle/l6as6ZTVt3kdx4wQVjiMicKxvPn4Ah2Mbjgw6E5/HVO3FxD+sfGI11HMCGGFw5g+BHypNObMIrehEJ92eR3nPVL9TWQ1H+HxzaUErf8qkwBWOIzpQ0PObNSXTF79fq+j9Civfj9l9W28ffCY11HMCGCFw5g+1OWeRFpHHRlt1V5H6VFOUwlZack8saXU6yhmBLDCYUwvjhxrpTVzCrn1+wf1lRI+7WLJaZNYtbOCts6A13HMMBfXwiEii0Vkn4gUicjNEeafIiJvikiHiNwUNq9ERHaKyDYR2RQyfYyIrBaRQvc+L57vwYxsT24tAyC3odDjJH27YsE0WjoDvLin0usoZpiLW+EQkSTgTuBSYD5wlYiEn2x+DLgRuK2HzVyoqmeo6sKQaTcDa1V1LrDWfW7MgFNVnthSSmZLGan+Zq/j9OmcmWOYmpvB41vKvI5ihrl47nGcAxSparGqdgIPA0tDF1DValXdCPhj2O5S4EH38YPApwYgqzHvseVwPSVHW8kdxAfFQ/l8wuVnTmV9YQ3Vje1exzHDWDwLx1Qg9MTyUndatBR4UUQ2i8i1IdMnqmoFgHs/IdLKInKtiGwSkU01NYkfbMcMff/+83uQoJ/RjcVeR4na5QumElRYsb3c6yhmGItn4Yh0LDGWk8zPV9UFOE1d14vIR2J5cVW9W1UXqurC8ePHx7KqMXQFgjTkzCan6RBJwVh2iL01Z3wW75uaw8odFV5HMcNYPAtHKTA95Pk0IOqfQapa7t5XA0/iNH0BVInIZAD3fvCeI2mGrLeKjxFIzmB0wwGvo8QkPz+fy06fwvYj9Rw51up1HDNMxbNwbATmisgsEUkFrgRWRLOiiGSKSHb3Y+DjwC539grgGvfxNcDTA5raGOCZHeX4Ap1kNx/2OkrMPnHaZACesb0OEydxKxyq2gXcALwAFACPqOpuEVkmIssARGSSiJQC3wZ+ICKlIpIDTATWi8h2YAPwrKo+7276VmCRiBQCi9znxgwYfyDI87sryW46hE+H3jUR08eM4ozpuTyzw45zmPiIa7fqqroKWBU2bXnI40qcJqxwjcD7e9jmUeCiAYxpzLu8ceAo9a1+ZjQWeR2l3y47fTI/e7aAg7UtzBpnQ8uagWVXjhsT5pnt5WSnJZPVPHR7m3301/8NOO/FmIFmhcOYEJ1dQV7YXcmiUyfi06DXcfottauFhTPy7DiHiQsrHMaEeL2olsb2Li47fbLXUY7bZadPZl9VE4VVTV5HMcOMFQ5jQqzcUU5OejIfOnHoX/uz5PTJiNjZVWbgWeEwxtXRFWD17iouOXUSqclD/7/GhOx0zp4xhhd2W6eHZmAN/f8dxgyQt4uP0dTRxYYn7h5UQ8Mej4+fOpG9lU18/IrPDZv3ZLxnhcMY15qCKjJSkshqGT69y15y6iQAGnNmeZzEDCd9Fg4ReV8ighjjJVVlzZ4qPjx33JC86K8n08eM4tQpOVY4zICKZo9juYhsEJHrRCQ33oGM8cKeikbKG9q5eP5Er6MMuMWnTqJ11CT8yaO8jmKGiT6vHFfVD4nIXODLwCYR2QA8oKqr457OmARZs6caVHngZ9+Ob3cKCZafn097Wh6c+G80Zs/0Oo4ZJqL6P6KqhSLyA2AT8DvgTBER4L9V9Yl4BjQm3vLz8ymafQWjNEhyYPgNgJTWUUdqR701V5kBE80xjtNF5A6cjgo/BuSr6jz38R1xzmdM3PmTM2nLmEBOU4nXUeJCgNGNxTRnTqGhdeiMLWIGr2iOcfwe2AK8X1WvV9Ut8M54GT+IZzhjEqExewYA2Y0l3gaJo5ymEpAk1u6t8jqKGQaiKRxPqOpfVLWte4KIfANAVf8St2TGJEhj9kxSOxpI66z3OkrcZLRVk+Jv5vlddjGgOX7RFI6rI0z74gDnMMYTLR1dtGROJaepJOJYx8OFADmNB3m1sIbWzi6v45ghrseD4yJyFfBZYJaIhI7clw0cjXcwYxLhtcJa1JdEdtMhr6PEXU7TQY6OPY3XCmvfuTDQmP7o7ayqN4AKYBxwe8j0JmBHPEMZkyhrCqpI6mons3X4N+FktlSSnZ7M2oIqKxzmuPRYOFT1EHAIOC9xcYxJnEBQeWlvNdnNhxGG7tgb0RKCXHjyBNYWVBMIKkm+4dw4Z+Kpx2McIrLevW8SkcaQW5OINEazcRFZLCL7RKRIRG6OMP8UEXlTRDpE5KaQ6dNFZJ2IFIjI7u6D8e68W0SkTES2ubclsb1lYxxbD9dxrKVzRDRTdbto3gSOtnSy7Ui911HMENZj4VDVD7n32aqaE3LLVtWcvjYsIknAncClwHzgKhGZH7bYMeBG4Law6V3Ad9zrRc4Frg9b9w5VPcO9rcKYGHT3Eru6oIqUJCF7CA8RG6sLTppAsk9YW2Cn5Zr+i+YCwDkikuY+vkBEboyyz6pzgCJVLVbVTuBhYGnoAqparaobAX/Y9IqQ60WacC4+nBrNGzImWmv2VHHu7LEkBTu9jpIwo0elcM6sMayxwmGOQzSn4z4OBETkROA+YBbw9yjWmwqE/pQrpR9//EVkJnAm8HbI5BtEZIeI3C8iebFu05jimmYO1LRw8bzh16lhXy6aN5H9Vc0cPtrqdRQzREVTOIKq2gVcDvxGVb8FRDMgc6QjbxpLOBHJwilc31TV7uMqdwFzgDNwzvq6vYd1rxWRTSKyqaamJpaXNSPA2oJqwGnzH2kudt+z7XWY/oqmcPjdazquAZ5xp6VEsV4pMD3k+TSgPNpgIpKCUzT+FtqRoqpWqWpAVYPAPThNYu+hqner6kJVXTh+/NAfP9oMrNUFVaS31TItb2R1NZ6fn88NX7ySuROyrHCYfoumcHwJ55Tcn6vqQRGZBfw1ivU2AnNFZJaIpAJXAiv6WAcAt+fd+4ACVf112LzQvZ3LgV3RbNOYbl1JaWwqOTZsOzWMxsXzJ7Lh4DEa2qzTQxO7PguHqu5R1RtV9SH3+UFVvTWK9bqAG4AXcA5uP6Kqu0VkmYgsAxCRSSJSCnwb+IGIlIpIDnA+8AXgYxFOu/2liOwUkR3AhcC3Yn/bZiRryjqBoDKiTsMNd/G8iXQFlVf2WzOuiV2f43GIyPnALcAMd3kBVFVn97Wue6rsqrBpy0MeV+I0YYVbT+RjJKjqF/p6XWN605g9k4k5aWS0j9w/mj+64RqSTr6aNXuq+OT7p3gdxwwx0QzkdB/Or/rNwPAZjNmMSB1dAZqzpnPZvInsfNPrNN4RlJymQ6zbl40/ECQlKZpWa2Mc0XxbGlT1Ofeai6Pdt7gnMyYO3io+RjAplUUj8DTccNlNJTS1d7Gx5JjXUcwQE03hWCcivxKR80RkQfct7smMiYM1e6qQoJ/z5oz1OornsptLSU32OeOtGxODaJqqPuDeLwyZpjhDxxozZKgqawqqyG4+QnpKktdxPOfTLtLqDvLXdXX88LJ5OCczGtO3PguHql6YiCDGxNvu8kYqGtqZNoLPpgqX3XSIpuwZFFY3c9LEbK/jmCEimr6qJorIfSLynPt8voh8Jf7RjBlYawqqEBnZp+GGy3E/C7sY0MQimmMcf8K5FqP7nL39wDfjlMeYuMjPz2f5yjc464Q8kgPtXscZNFK6Wkhvq3mnCxZjohFN4Rinqo+AM9KNe2GfnZZrhhR/cibtGeO5yM6meo+cphK2HK6jtrnD6yhmiIimcLSIyFjcDgpF5FygIa6pjBlgjdkzAFg0f+R1atiXnKZDqMK6vbbXYaITTeH4Nk4fU3NE5HXgz8B/xDWVMQOsMXsmqR0NzBmf5XWUQSe9vZZJOel2nMNELZqzqraIyEeBk3G6AdmnqtYzmhkS8vPzCfhSaDn5i4w9tstOOY1AcLqXf3JrGe3+gJ2qbPrUY+EQkSt6mHWSiBDa1bkxg1lz5jTUl0T2CO4Nty/rH1lO64wlvFV8lAtOtuY807ve9jjy3fsJwAeBl9znFwIvA1Y4zJDQmDOTpK52MlsrvY4yaGW1lCFBP2sKqqxwmD71WDhU9UsAIvIMMF9VK9znk4E7ExPPmOOjCE1ZJ5DdfBiJbQDKEcWnAbKbS3mpIBtdqtakZ3oVzcHxmd1Fw1UFnBSnPMYMqNZREwkkZ4zoQZuild1UQnlDO3sqGvte2Ixo0fRV9bKIvAA8hHNK7pXAurimMmaANGbPRIIBspqPeB1l0MtuOowIrNlTzalTRnsdxwxi0YwAeAOwHHg/cAZwt6ra6bhmSGjMnklmazlJQTsRsC8pgTbOmJ7L2r12Wq7pXTR7HKjqk8CTcc5izIA6UNNMZ1ou447t9DrKkHHxvIn86oV9VDW2MzEn3es4ZpCyYb/MsLVmj/PL2To1jN5Td/4EwPquMr2Ka+EQkcUisk9EikTk5gjzTxGRN0WkQ0RuimZdERkjIqtFpNC9z4vnezBD15qCKtLbakn1N3sdZchI66gjpbORtXYVuelFNN2qXyYiMRcYEUnCOW33UmA+cJWIzA9b7BhwI3BbDOveDKxV1bnAWve5Me9ytLmDzYfq7GyqGAlO31Xri2pp67S+TE1k0RSEK4FCEfmliMyLYdvnAEWqWqyqncDDwNLQBdxxzDcC4Ucue1t3KfCg+/hB4FMxZDIjxLp9NQQVKxz9kNNUQkdXkPVFtV5HMYNUNGdVfR44EzgAPOA2LV0rIn0NFzYVCD0HstSdFo3e1p3YfV2Je2+XuZr3WLOniok5aaS32x+/WI1qrSA7Ldmaq0yPomqCUtVG4HGcX/6TgcuBLSLS22m5kS49jfbS3eNZ19mAU9w2icimmpqaWFY1Q1y7P8CrhTVcPG9ixC+S6Z1Pg3zk5PGsKagmGLSr7c17RXOM45Mi8iROX1UpwDmqeinOdR039bJqKTA95Pk0oDzKXL2tW+V2e9Ld/UnE0z9U9W5VXaiqC8ePHx/ly5rh4M3io7R2Brh4vg3a1F8Xz5tAbXMHO8ps6B3zXtHscfwLcIeqnq6qv1LVagBVbQW+3Mt6G4G5IjJLRFJxjpWsiDJXb+uuAK5xH18DPB3lNs0IsWZPFaNSkzhv9livowxZF5w0AZ/885RmY0JFUzgqVPXV0Aki8gsAVV3b00ruELM34IxXXgA8oqq7RWSZiCxztzNJREpxBov6gYiUikhOT+u6m74VWCQihcAi97kxAASDypqCKpJrC/nXKz7ldZwhKy8zlYzmchvcyUQUzZXji4DvhU27NMK091DVVcCqsGnLQx5X4jRDRbWuO/0ocFGfqc2ItLOsgarGDqY1lngdZcjLbjrE3soplNa1Mi1vlNdxzCDS4x6HiHxdRHYCp4jIjpDbQWBH4iIaE701BVX4BLKbD3sdZcjrPpXZriI34Xprqvo7zmBOT7v33bez3FN0jRl0Vu+pYuHMMSQH2r2OMuSldTYwe1ymNVeZ9+itcKiqlgDXA00hN0RkTPyjGRObI8da2VvZxMftbKoBc9G8CbxVfJSmdutd2PxTX3scAJuBTe795pDnxgwqL7pnAC2ywjFgLp43EX9Aea3QLqQ0/9Tb0LGXufezEhfHmP5bvaeSuROymDE20+sow8ZZM/JI6mpnTUEVS06b7HUcM0j0WDhEZEFvK6rqloGPY0z/1Ld2srGkjjFVm8nPv63vFUxUkpN8ZDcfZt3ebAJBJcln1+Kb3k/Hvb2XeQp8bICzGNNvL+2tJhBU69QwDrKbDnGk9SS2HK7j7Jl2eNP03lR1YSKDGHM8brnnCZJHTSKjzU4dHWjZzUdI9glrCqqscBig96aqj6nqSyJyRaT5qvpE/GIZE72OrgDNWdPJbSiyTg3jICnYybmzx7JmTxX/dWksIyuY4aq3s6o+6t7nR7hdFudcxkTtjQNHCSalkm3NVHFz0bwJHKhp4WBti9dRzCDQW1PVj9z7LyUujjGxW72nCl/AT1ZLmddRhq2L503kxyv3sLagiq9+eLbXcYzHoulWfayI/E5EtojIZhH5rYhYt6NmUAgGlTV7qshqOYJPbajTeJk+ZhRp7Ue545Ee+zU1I0g0veM+DNQAn8bpYr0G+Ec8QxkTrR1lDVQ3dZBjnRrG3ejGg7SMmkxNU4fXUYzHoikcY1T1p6p60L39DMiNcy5jovLC7kqSfEJ28yGvowx7OY0HQYTVNkbHiBdN4VgnIleKiM+9fQZ4Nt7BjOmLqnLfC5tJbzxCcsB+BcdbesdRUjsaeH53pddRjMd6Ox23CedCP8EZaOmv7iwf0Az8KO7pjOnFvqomOtNyGXfUevlPBAFymg7yRlEuDW1+RmekeB3JeKTHPQ5VzVbVHPfep6rJ7s2nqjmJDGlMJM/trARVcpoOeh1lxBjdWExXUFlrXa2PaNE0VSEieSJyjoh8pPsW72DG9OWF3ZWMaq0kpavN6ygjRkZbNZNy0nl+lzVXjWR9Dh0rIl8FvoEzxOs24FzgTayvKuOhg7Ut7K1sYnJjsddRRhQBFr9vEg9tOExLRxeZadGMPm2Gm2j2OL4BnA0ccvuvOhPnlNw+ichiEdknIkUicnOE+eJeI1LkDku7wJ1+sohsC7k1isg33Xm3iEhZyLwl0b5ZM3xc+d1bAayZygOXnDqJjq4gr+yP6s+AGYaiKRztqtoOICJpqroXOLmvlUQkCbgTuBSYD1wlIvPDFrsUmOvergXuAlDVfap6hqqeAZwFtAJPhqx3R/d8VV0VxXsww0xjzmwyWqtJ9Td7HWXEOWfWGMZmpvKcNVeNWNEUjlIRyQWeAlaLyNNAeRTrnQMUqWqxqnbiXEi4NGyZpcCf1fEWkCsi4aPFXAQcUFU7Ud8AUFrXSlvGBEY3WTOVF5J8wqL5E3mpoIp2v12tPxL1WThU9XJVrVfVW4AfAvcBn4pi21OBIyHPS91psS5zJfBQ2LQb3Kat+0UkL9KLi8i1IrJJRDbV1Ngu9XDSfWA2x45veCI/P5+3HltOS2eA14tsSNmRKNqzqhaIyI3A6UCpuwfR52oRpmksy4hIKvBJ4NGQ+XcBc4AzgAp6GHBKVe9W1YWqunD8+PFRxDVDxfO7KklvryWts9HrKCNWZksZ2enJdnbVCBVNJ4f/AzwIjAXGAQ+IyA+i2HYpMD3k+TTe28TV1zKXAltU9Z2TxlW1SlUDqhoE7sFpEjMjRHVjO5sP1zndXxjP+DRIUuUenni7CH8g6HUck2DR7HFcBZytqj9yu1o/F/hcFOttBOaKyCx3z+FKYEXYMiuAq92zq84FGlS1Iuy139VMFXYM5HJgVxRZzDCx9D9+jKpzIZrx1uiGIgLJ6awvtOaqkSaawlECpIc8TwMO9LWSqnYBNwAvAAXAI6q6W0SWicgyd7FVQDFQhLP3cF33+iIyClgEhI80+EsR2SkiO4ALgW9F8R7MMNEw+kROmZRNeked11FGvKyWUnyBDlZuj+ZcGTOc9NZX1f/hHG/oAHaLyGr3+SJgfTQbd0+VXRU2bXnIYwWu72HdVpzmsfDpX4jmtc3wU1rXSuuoSeS/fwrP27AQnvNpkNGNxby4J5N2f4D0lCSvI5kE6e2yz03u/WbefQ3Fy3FLY0wvnt3htGLmnz6F5z3OYhy5DUUczJvHy/uqWfy+8DPpzXDV29CxD3Y/do9RnOQ+3aeq/ngHMybcyh3lZLRWccLYUV5HMa7MlnLGZaWycnuFFY4RJJqzqi4ACnGuAv8DsN86OTSJVlzTzK6yRnIb+zy8ZhJIUJacNpm1e6to7ujyOo5JkGgOjt8OfFxVP6qqHwEuAe6Ibyxj3m3l9gpEYHSDFY7BZv3ff0u7P8ilX/ku+fn5XscxCRBN4UhR1X3dT1R1P2AjuJiEUVVWbC/j7JljSOlq8TqOCTOqtZIUfzP1OSd6HcUkSDSFY7OI3CciF7i3e3AOmBuTEHsrmzhQ00L++6d4HcVEIDjXdDRnTSfgS/U6jkmAaArHMmA3cCNOF+t73GnGJMTK7eUk+YQl75vkdRTTg9ENB1BfEg05s72OYhKg11FYRMQHbFbV9wG/TkwkY/4pGFSe3lZORsMhxmaleR3H9CCjvYbUjgbqR8/1OopJgF73ONz+oLaLyAkJymPMu7x98Bhl9W3kNez3OorphQC5DftoyZpKaV2r13FMnEXTVDUZ58rxtSKyovsW72DGADy+pZSstGRyGku8jmL6kFdfCMBTW8s8TmLiLZoBg38c9xTGRNDa2cVzOyv4xOmTKdhi1wgMdqn+JjJbynl8SybXX3giIpFGTTDDQY97HCKS7o7z/a/AKcDrqvpK9y1RAc3I9cLuSlo6A3x6wTSvo5go5dbv42BtC1uP1HsdxcRRb01VDwILgZ0442JEHDDJmHh5YksZ0/IyOHvmGK+jmCiNbiwmPcXH45tLvY5i4qi3wjFfVT+vqn8E/gX4cIIyGUNFQxvri2q54syp+HzW5DFUJAX9XHLqJFZuL7fxyIex3grHOx0ZumNrGJMwT20tRxWe/+NPvI5iYrT9qbtpbO/ipb3VXkcxcdLbwfH3i0j3oM4CZLjPBWcojZy4pzMj0mX5+RTO+Qyjgp02rvgQlNVSRrK/hcc3l7LkNOsxdzjqrVt1G5XFeKI1YxId6WOYWvay11FMPwhKXv1+Xt6fRWVDO5NGp/e9khlSormOw5iEOpY3D1+gk9GNRV5HMf2UV19AIKg8uumI11FMHFjhMINKQ6ufhtFzyG0oJCloh9aGqrTORjKbS/nds5sJBNXrOGaAxbVwiMhiEdknIkUicnOE+SIiv3Pn7xCRBSHzSkRkp4hsE5FNIdPHiMhqESl07/Pi+R5MYj21rQz1JTOmrsDrKOY4jakrwJ+azWuFNV5HMQMsboVDRJJwRg28FJgPXCUi88MWuxSY696uBe4Km3+hqp6hqgtDpt0MrFXVucBa97kZBlSVhzYcJr2thoz2Wq/jmOOU03SQpK42Htpw2OsoZoDFc4/jHKBIVYtVtRN4GFgatsxS4M/qeAvIFZG+TsNYinNxIu79pwYws/HQtiP17K1sYkzdHq+jmAHg0yB59ftYU1BNdWO713HMAIpn4ZgKhB4ZK3WnRbuMAi+KyGYRuTZkmYmqWgHg3k+I9OIicq2IbBKRTTU1tqs8FDy84QijUpPIbbCD4sPFmDr3ILldST6sxLNwRLrcN/woWW/LnK+qC3Cas64XkY/E8uKqereqLlTVhePHj49lVeOBhlY/K7aX88n3TyEp6O97BTMkpHU28ME5Y/nNyo12kHwYiWfhKAWmhzyfBpRHu4yqdt9XA0/iNH0BVHU3Z7n3dnnqMPCPTYdp8we4+ryZXkcxA+xzH5iBPzWHtQVVXkcxAySehWMjMFdEZolIKnAlED6OxwrgavfsqnOBBlWtEJFMEckGEJFM4OPArpB1rnEfXwM8Hcf3YBIgEFQefOMQmS3lfO9rn/M6jhlgl5w6kZTOJh54vcTrKGaAxK1wuP1b3QC8ABQAj6jqbhFZJiLdY5avAoqBIuAe4Dp3+kRgvYhsBzYAz6rq8+68W4FFIlIILHKfmyHsws/fSFl9G2OP7vQ6iomD5CQfY+t282bxUQoqrAuZ4SCagZz6TVVX4RSH0GnLQx4rcH2E9YqB9/ewzaPARQOb1Hjp6JjTSOlsIqepxOsoJk7y6gqon/ZBPnvL3Wz9401exzHHya4cN57aW9lIS9ZUxh7bjbzn3AkzXCQHOrhiwTTqR8/lWEun13HMcbLCYTz1p9dLkKCfvHq7Uny4+9IHZ6K+ZLsgcBiwwmE8U93YzhNbysir309yoMPrOCbO5k7MJqv5CH9+s4SOLhvkaSizwmE8s+Tbt9HZ1cW42u1eRzEJMq52O1WNHTy1tczrKOY4WOEwnmhs93Msbz6jG4tJ89uZNiNFVksp75uaww///hqX5X/S6zimn6xwGE/89a1DBJPSGF+71esoJoEEuO6CE+lMy6UhZ7bXcUw/WeEwCdfuD3D/+hKymg6T0X7U6zgmwS45dRJpHXXUjDsT54x8M9RY4TAJ9+imI9Q2dzC+dpvXUYwHknzCuNpttGeM45X91gHpUGSFwyRUuz/AnesOcNaMPDJbw7suMyNFbkMhKZ1N/N9LRbbXMQRZ4TAJ9fe3D1PZ2M53Fp0UsWtkMzL4NMj42q1sPlTHy7bXMeRY4TAJ09YZ4A8vH+Dc2WP44InjvI5jPJZXv5dpeRnc/uI+2+sYYqxwmIT585sl1DZ3UL32fvLz872OYzzm0yDsfJZdZY08v6vS6zgmBlY4TEI0d3Txy5VbyWo+Qmar/ZEwjtyGQtLaj3H76v020NMQYoXDJMRdLxcRSM5gYtUGr6OYQURQJtZsoqi6mSftavIhwwqHibuy+jbufe0gufX7GdVuB0LNu+U0FvP+aaO57YV9tHZ2eR3HRMEKh4m7Xz2/F4CJ1ba3Yd5LgP/Jn09lYzvLXyn2Oo6JghUOE1fbjtTz1LZyvvrhWaT6m72OYwapW264htENhfzxlQOU1bd5Hcf0wQqHiZtgUPnxyt0kd7Xy9QtO9DqOGeQmVb0NwC+e2+txEtMXKxwmbh7eeISth+uZVPkWWWlxHaXYDAOp/ma+9pHZrNhezpsHrA+zwSyuhUNEFovIPhEpEpGbI8wXEfmdO3+HiCxwp08XkXUiUiAiu0XkGyHr3CIiZSKyzb0tied7MP1T09TBrc8VcO7sMeQ27Pc6jhki1v7+ZlI7G/j+kztp99tgT4NV3AqHiCQBdwKXAvOBq0RkfthilwJz3du1wF3u9C7gO6o6DzgXuD5s3TtU9Qz3tipe78H03/9bVUCbP8DPPnWadS1ioubTLqaUv0ZxbQt/WFfkdRzTg3jucZwDFKlqsap2Ag8DS8OWWQr8WR1vAbkiMllVK1R1C4CqNgEFwNQ4ZjUD6NX9NTy5tYxlH53DiROyvI5jhpjsllJy6/fxf2v2sb+qyes4JoJ4Fo6pwJGQ56W8949/n8uIyEzgTODtkMk3uE1b94tIXqQXF5FrRWSTiGyqqbFrBxKlodXPfz62gxMnZPHS77/ndRwzRE2ufBNfsJP/fGwHXYGg13FMmHgWjkgtFOF9CvS6jIhkAY8D31TV7vFF7wLmAGcAFcDtkV5cVe9W1YWqunD8+PExRjf99aMVu6ht7uCOz5yBT62N2vRPcqCdKRXr2Xaknv97yZqsBpt4Fo5SYHrI82lA+AAMPS4jIik4ReNvqvpE9wKqWqWqAVUNAvfgNImZQWDVzgqe2lbOmMoNnDZttNdxzBCX23iAK86cyu/XFbH5UJ3XcUyIeBaOjcBcEZklIqnAlcCKsGVWAFe7Z1edCzSoaoWICHAfUKCqvw5dQUQmhzy9HNgVv7dgonXkWCs3P76DjNZqJtTYOOJmYOz9+0/xtTdy1W+fY8nSK7yOY1xxKxyq2gXcALyAc3D7EVXdLSLLRGSZu9gqoBgowtl7uM6dfj7wBeBjEU67/aWI7BSRHcCFwLfi9R5MdDq6Anz8x/+guaWFE0pXI1ibtBkYScFOppetpTMlm9IpF9i4HYNEXK/Kck+VXRU2bXnIYwWuj7DeeiIf/0BVvzDAMc1x+tkzBbRlTGDG4edJ9dtZMGZgZbZWMqnqbSonncc9rxVz7UfmeB1pxLMrx81xeXTTEf7y1iHG1W4jp6nE6zhmmBp3dDs5DQe49bm9vFFU63WcEc8Kh+m3Nw7U8t1HncGZJtk4GyaOBJhW/jKzx2dx3d+3cKDGOsz0khUO0y9F1c0s+8tm0joaOOGIHdcw8ZcU9HPfNQtpamzkiw9soLqp3etII5YVDhOz8vo2rrl/AylJPmYefo6kYKfXkcwIMWNsJjMPraK2qZOP/uAhliy93OtII5IVDhOTxZf/Gxf8+Eka2/z86Uvn2MFwk3Cj2mv4/WfPpC19HIdOWEJLh40amGhWOEzUapo6ODjjMrqSRzG+4FH+++uf9zqSGaEumjeR6aVraRk1iXNuesD2PBLMCoeJypFjrfzr8jfoTM1mxuHnyGyr8jqSGeFyGw+8UzxKTlhCQ5vf60gjhhUO06d9lU1c+PNnOFJ1jFklz5DVGt5zjDHe6C4ebRkT+Ze73rBhZxPECofp1eo9VXz6rjcAmF3ytO1pmEEnt/EAMw8/y4GKo3z0pyvZWdrgdaRhzwqHiSgYVH63tpB///MmZo3LZE7xE6R3WEdzZnDKailnzsGnEA3yyd+9zPlX/6d1TxJHVjjMe1Q2tHP1/Rv49er95NbvI7j6NlK7WryOZUyv0jvqOLH4cTJbKyib8lG+8+h2mtrtuEc8WOEw71BVVm4v55LfvMrr+yv4f5efxrSydTauhhkykgPtzDy0ignVm3hqaxmLf/MarxXaQG4DzQqHAaC4pplrHtjIfzy0lVnjMpl74DE++4ETbLxwM+QIysSaTcw88CS1VeV84b4NfPfR7Xal+QCywjHC1bV08r/PFbD4N6+x9VAdkyvW89iy80jrtAOMZmjLbKvixAOPseyjc3hqWxkX/upl7lxXRLvf9qCPlxWOEaqupZPbX9zHh37xEn98+QCjjhaw9qaPMu7YLpKT7GthhgefBnj9D99l1t6/k1RbxK9e2MeHfrGOP75ygGa74rzf4joehxl89pQ38uAbJTy1rYyOriCfOH0y+5/4DekddUzI/rbX8YyJi7TOBmYceYGWUZOY+vGv8r/P7eUPLx/gynOmc+XZJzBrXKbXEYcUKxwjQFVjOyu3l7NyeznbSxuQoJ+8+kJOOLaTw7vrSPc6oDEJktlayV++8gEu+sxXmJt/Hfe+dpA/vlLMB+eM5fIzp7Jo/kRyR6V6HXPQs8IxDAWDyu7yRl4trOHlfdVsOlSHKqS31TC5oZC8ur3Wo60Z0Ua1VbP8C2ex+PJ/I//Gn/PwxsN897EdJPmE82aPZdH8iZx/4ljmjM9CxE4RCWeFYxhoavezs7SBrUfq2Xq4nnU7SwgkZwBOsRjfdIjchkI74G1MmJSuVq6/8ESe+/W3SE0fT2POLDa0zWa9O8pgsr+FzJZyvvH5pZw+bTTzp+QwKtX+bMb1ExCRxcBvgSTgXlW9NWy+uPOXAK3AF1V1S2/risgY4B/ATKAE+IyqDvtLmls7u6hoaKeyoZ2Soy38+t6/05GaS0daHv7U7HeWmz0uk+zmw2Q1l5LVXEpKwPruMaYvgtNd+6j2GiZWb8Cfkk1z5lSaM6fSkjmFnzyzBwCfwMxxmcwel8Wc8ZnMGufcpuRmMCEnjbTkJG/fSILErXCISBJwJ7AIKAU2isgKVd0TstilwFz39gHgLuADfax7M7BWVW8VkZvd59+L1/s4XqpKZyBIuz9Ihz9Auz9Ie1eA9u7H/gCtnQEa2/w0tPmpb+t07lud5xt27sOfkkUwKe1d25W8eaR11DOqtZL0ugJ+cfP1/OJ71/HcU4+Rn3+7R+/WmKFPgFR/E2Pq9zKmfi8A/uRRtKWPpy1jPLUNY0n2ncarhTV0dr175MsxmalMyE5j0uh08kalMjojhZz0ZHIyUshJT3Hvk0lPTSIjJYn0lO57H+kpSST7hCSfDPrmsXjucZwDFKlqMYCIPAwsBUILx1Lgz+p0KvOWiOSKyGScvYme1l0KXOCu/yDwMnEqHMtfOcCzOyoIBJWgKoGgElAl+M4970x/Z35QCSrvLOsPBImpyxwNkhTo4IRJ48jJSCGts5GslnJSulpI9reQ4m8m1d9Eir/5XRfnXXDyBG4PdAz0R2CMwWnSSmk+RE7zIWfCkRc4CeGuBx/mC9d9h+u+89/cdte9LMr/F1asfoVDKZl0JaUT9KUSSEqDGAtBSpKQ7PORnCSkJPlI9jn3Ph/4RPCJIOI8FnjXc58PBMEnICL88LL5nDUjb0A/j3gWjqnAkZDnpTh7FX0tM7WPdSeqagWAqlaIyIRILy4i1wLXuk+bRWRfL1nHAbW9zE+44siTe8zZ/QvleH6phG+jn9saB9QO0LYGMlf4Nt6VcyByDcQ2wrYV0/fSw8+815zxzBWDdzIOZJ4Txjqn8b7yj+UAbHzgJ/3elmvA/xY9fcNxrT4j0sR4Fo5I/xrhv717WiaadXulqncDd0ezrIhsUtWFsWzfC0Mh51DICEMj51DICEMj51DICEMnZzwvES4Fpoc8nwaEjwDU0zK9rVvlNmfh3lcPYGZjjDF9iGfh2AjMFZFZIpIKXAmsCFtmBXC1OM4FGtxmqN7WXQFc4z6+Bng6ju/BGGNMmLg1Valql4jcALyAc0rt/aq6W0SWufOXA6twTsUtwjkd90u9retu+lbgERH5CnAY+NcBiBtVk9YgMBRyDoWMMDRyDoWMMDRyDoWMMERyio2SZYwxJhbWDaoxxpiYWOEwxhgTk2FdOETkfhGpFpFdIdN+JSJ7RWSHiDwpIrk9rFsiIjtFZJuIbPIg5y0iUua+/jYRWdLDuotFZJ+IFLlX0icy4z9C8pWIyLYe1k3IZyki00VknYgUiMhuEfmGO32MiKwWkUL3PuLVUAn8LHvKOWi+m71kHGzfy55yDprvpoiki8gGEdnuZvyxO31QfS9joqrD9gZ8BFgA7AqZ9nEg2X38C+AXPaxbAozzMOctwE19rJcEHABmA6nAdmB+ojKGzb8d+B8vP0tgMrDAfZwN7AfmA78Ebnan3xzp3zzBn2VPOQfNd7OXjIPtexkx52D6buJcl5blPk4B3gbOHWzfy1huw3qPQ1VfBY6FTXtRVbuH/noL5xoRT0XKGaV3unVR1U6gu2uWAddbRhER4DPAQ/F47WipaoW6nWSqahNQgNMLwVKc7mlw7z8VYfVEfpYRcw6m72Yvn2U0PP8su+cPhu+mOprdpynuTRlk38tYDOvCEYUvA8/1ME+BF0Vkszjdl3jhBrfZ4v4edmN76rIl0T4MVKlqYQ/zE/5ZishM4EycX3fv6qYGiNRNjSefZVjOUIPmuxkh46D8XvbwWQ6K76aIJLnNZdXAalUd1N/LvozYwiEi3we6gL/1sMj5qroApwff60XkIwkL57gLmAOcAVTg7G6HO+6uWQbIVfT+iy6hn6WIZAGPA99U1cZoV4swLa6fZU85B9N3M0LGQfm97OXffFB8N1U1oKpn4OxFniMi74ty1cHyf/xdRmThEJFrgMuAz6nbkBhOVcvd+2rgSZxdxoRR1Sr3yxYE7unh9aPp1iWuRCQZuAJnjJSIEvlZikgKzh+Qv6nqE+7kaLqpSehn2UPOQfXdjJRxMH4ve/ksB9V3032depwevRczCL+X0RpxhUOcAaK+B3xSVVt7WCZTRLK7H+MctNwVadl46f5CuS7v4fWj6dYl3i4G9qpqaaSZifws3fbs+4ACVf11yKxouqlJ2GfZU87B9N3sJeOg+l728m8Og+S7KSLjxT1DTkQyunMxyL6XMfH66Hw8bzi7qBWAH6dyfwWne5MjwDb3ttxddgqwyn08G+fshe3AbuD7HuT8C7AT2IHzRZkcntN9vgTnTJID8cwZKaM7/U/AsrBlPfksgQ/h7MbvCPn3XQKMBdYChe79GI8/y55yDprvZi8ZB9v3MmLOwfTdBE4HtroZd+Ge4TXYvpex3KzLEWOMMTEZcU1Vxhhjjo8VDmOMMTGxwmGMMSYmVjiMMcbExAqHMcaYmFjhMMYYExMrHMYYz4jIPBFZLiKPicjXvc5jomOFwySEiKiI/CXkebKI1IjIM17mGmgi8oZ7nysi18Vr+73MnykhY6aETO81j7temzsuRcRtuMupiNwe8vwmEbklbJk/isj5fb4ZQFULVHUZTg+2C931M9wcnSIyLprtmMSywmESpQV4n9vlAsAioCzRIcQRt++9qn7QfZgLDHjhCNl+rHLpO88BdTri600HcEUff9A/gNMtfFRE5JPAepyrp1HVNjeH530ymciscJhEeg74hPv4Xb2WishTbtfWu7u7t3b7EnpWnJHTdonIv7nTv+0+3yUi3wzZxkxxRoK7x93Oi+6v1+7pfwC2APeJO1Kcu97PReTG0KDhv7q7f1n39Bohy3WPu3ArMMf95fyrsG3/Z/fricgdIvKS+/giEfmr+/jz4owat839BZ8Utn1E5IfijBi4WkQeEpGb3FlJEfL1mKc3IjJbRLaKyNnupC7gbuBbPSw/D6d7jOlutnvdf6e/icjFIvK6OCPevdOZoKqucAvi56LNZTzmdZ8ndhsZN6AZp8+ex4B0nD6FLgCeced399OTgdOfz1jg08A9IdsYDZyF01dSJpCF08fQme78mTh/2M5wnz8CfN6dHgTODVlui/vYh9MH0NiwvDN594iMN+GMfhfxNULfZ6T1w7Z9LvCo+/g1YAPO4D4/Ar4GzANWAinuMn8Arg7b/kL3M8zAGfmu0M3Y22cQMU943u7HwMk4fSydEfbvmIMzct7o7s8lZP63ccYS6c5xmvsZbwbux+kmfCnwlLv8BcDvgD8C14dlKiFBo3DaLbZbMsYkiKruEGewnauAVWGzbxSRy93H04G5OAXiNhH5BU6BeU1EPgQ8qaotACLyBM5gPVvddQ+q6jb38WacP2DrgUOq+pabo0REjorImcBEYKuqHo3hrUR6jVhsBs4Sp2fWDpy9oIXu+7gRuAinQG4UEXCKQ3iX2x8CnlbVNgARWdlHvvUxZhyP01vrp1V1d+gMVW0UkT+7WdvC1rsE+BLOMKcHVXWnm283sFZVVUR2uplQ1Zdxuhk3Q4gVDpNoK4DbcH5pjgUQkQtwupo+T1VbReRlIF1V94vIWTi9g/6viLwI9DUwU0fI4wDOH11wjrGEuhf4IjAJ55dwuC7e3ZSbHsVrREVV/SJSgvMH9g2cXlMvxBkgqQCncDyoqv/Vy2YiDfAzIPlcDTg99Z6Ps1cX7jc4Be+BdwKJjAJyVbXc/YEQmiMY8jyI/e0Z0uwYh0m0+4GfdP8SdY0G6tyicQpOUw4iMgVoVdW/4hSbBcCrwKdEZJQ4YyhcjtPcE6sncQbTORt4IcL8KmCCiIwVkTScwZVi0YTThNSTV3GaeV7Fyb8M2KaqinOQ+F9EZAKAiIwRkRlh668H8kUkXZzR7z5B7/rKE64TZwzsq0Xks+EzVfUYTjPYV0ImXwisi+E1zBBlhcMklKqWqupvwyY/DySLyA7gp/zzjJzTgA3ijNX8feBnqroFZ5yFDThjS9+rqluJkap24vyRe0RVAxHm+4GfuK/xDM7AO7Fs/yjwuntgONLB6NeAycCbqloFtLvTUNU9wA9wxsLeAax2lw3d/kacvbftwBPAJpy9hP7mibROC07B/JaILI2wyO1A6NlVl+L8W5phzsbjMCOSe0ruFuBfVbXQ6zz9ISJZqtrsNhG9ClzrFtb+bGsmznGkaMfCjrSNLcAH3KJ73NzmvIWqWjsQ2zMDx/Y4zIgjIvNxRttbO1SLhutud29sC/B4f4uGKwCMdrfXL6q6YCCKhnsK9TacM82Cx7s9M/Bsj8MYY0xMbI/DGGNMTKxwGGOMiYkVDmOMMTGxwmGMMSYmVjiMMcbExAqHMcaYmFjhMMYYExMrHMYYY2JihcMYY0xM/j9pS7qqtZpg6QAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#https://realpython.com/python-histograms/\n",
    "n, bins, patches = plt.hist(x=r, bins='auto', color='black', alpha=0.7, rwidth=0.85, density=True)\n",
    "\n",
    "plt.plot(np.sort(r), compressive_resistance_distribution)\n",
    "plt.xlabel('Masonry unit weight [kN/m$^3$]')\n",
    "plt.ylabel('Probability density')\n",
    "# plt.xlabel(x_label_pdf)\n",
    "# plt.ylabel(y_label_pdf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 114,
   "metadata": {},
   "outputs": [],
   "source": [
    "realisations_vector = np.transpose(np.sort(r))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 115,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "100000"
      ]
     },
     "execution_count": 115,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.sort(r).shape[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 118,
   "metadata": {},
   "outputs": [],
   "source": [
    "realisations_series = pd.Series(realisations_vector)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 119,
   "metadata": {},
   "outputs": [],
   "source": [
    "realisations_header = 'masonry_unit_weight_realisations'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 122,
   "metadata": {},
   "outputs": [],
   "source": [
    "realisations_vector_file_name = 'masonry_unit_weight_realisations.csv'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 129,
   "metadata": {},
   "outputs": [],
   "source": [
    "realisations_series.to_csv(realisations_vector_file_name, header = [realisations_header], index = False)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Moments calculation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {},
   "outputs": [],
   "source": [
    "second_moment_list = []"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {},
   "outputs": [],
   "source": [
    "# for i in range(simulations_number):\n",
    "#     print(i)\n",
    "#     print(realisations_list[i])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "metadata": {},
   "outputs": [],
   "source": [
    "def second_moment_calculation(realisations):\n",
    "    for i in range(simulations_number):\n",
    "        sample_second_moment=stats.moment(realisations[i], moment=2)\n",
    "\n",
    "        second_moment_list.append(sample_second_moment)\n",
    "\n",
    "        #print(second_moment_list)    #return second_moment_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "metadata": {},
   "outputs": [],
   "source": [
    "second_moment_calculation(realisations_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "metadata": {},
   "outputs": [],
   "source": [
    "sample_mean_cov_list = []\n",
    "r_array = np.array(r)\n",
    "sample_mean = np.mean(r_array)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 99,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sample_mean_cov_calculation(second_moments):\n",
    "\n",
    "    for i in range(simulations_number):\n",
    "        sample_mean_cov = second_moments[i]/sample_mean\n",
    "        \n",
    "        sample_mean_cov_list.append(sample_mean_cov)\n",
    "\n",
    "    #return sample_mean_cov_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 100,
   "metadata": {},
   "outputs": [],
   "source": [
    "sample_mean_cov_calculation(second_moment_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "metadata": {},
   "outputs": [
    {
     "ename": "ValueError",
     "evalue": "x and y must have same first dimension, but have shapes (144,) and (100000,)",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[0;32m/var/folders/x7/6bdwhg5j1lg_m9s3kq825p380000gq/T/ipykernel_31811/891838005.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mplot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mn\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msample_mean_cov_list\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcolor\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'black'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      2\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mxticks\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0msimulations_number\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mxlabel\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"Number of simulations\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mylabel\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"Coefficient of variation of sample mean\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/opt/anaconda3/lib/python3.9/site-packages/matplotlib/pyplot.py\u001b[0m in \u001b[0;36mplot\u001b[0;34m(scalex, scaley, data, *args, **kwargs)\u001b[0m\n\u001b[1;32m   3017\u001b[0m \u001b[0;34m@\u001b[0m\u001b[0m_copy_docstring_and_deprecators\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mAxes\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mplot\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   3018\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0mplot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mscalex\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mTrue\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mscaley\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mTrue\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdata\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mNone\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 3019\u001b[0;31m     return gca().plot(\n\u001b[0m\u001b[1;32m   3020\u001b[0m         \u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mscalex\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mscalex\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mscaley\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mscaley\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   3021\u001b[0m         **({\"data\": data} if data is not None else {}), **kwargs)\n",
      "\u001b[0;32m/opt/anaconda3/lib/python3.9/site-packages/matplotlib/axes/_axes.py\u001b[0m in \u001b[0;36mplot\u001b[0;34m(self, scalex, scaley, data, *args, **kwargs)\u001b[0m\n\u001b[1;32m   1603\u001b[0m         \"\"\"\n\u001b[1;32m   1604\u001b[0m         \u001b[0mkwargs\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcbook\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mnormalize_kwargs\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmlines\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mLine2D\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1605\u001b[0;31m         \u001b[0mlines\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_get_lines\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdata\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mdata\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1606\u001b[0m         \u001b[0;32mfor\u001b[0m \u001b[0mline\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mlines\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1607\u001b[0m             \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0madd_line\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mline\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/opt/anaconda3/lib/python3.9/site-packages/matplotlib/axes/_base.py\u001b[0m in \u001b[0;36m__call__\u001b[0;34m(self, data, *args, **kwargs)\u001b[0m\n\u001b[1;32m    313\u001b[0m                 \u001b[0mthis\u001b[0m \u001b[0;34m+=\u001b[0m \u001b[0margs\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    314\u001b[0m                 \u001b[0margs\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0margs\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 315\u001b[0;31m             \u001b[0;32myield\u001b[0m \u001b[0;32mfrom\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_plot_args\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mthis\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    316\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    317\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0mget_next_color\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/opt/anaconda3/lib/python3.9/site-packages/matplotlib/axes/_base.py\u001b[0m in \u001b[0;36m_plot_args\u001b[0;34m(self, tup, kwargs, return_kwargs)\u001b[0m\n\u001b[1;32m    499\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    500\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mx\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m!=\u001b[0m \u001b[0my\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 501\u001b[0;31m             raise ValueError(f\"x and y must have same first dimension, but \"\n\u001b[0m\u001b[1;32m    502\u001b[0m                              f\"have shapes {x.shape} and {y.shape}\")\n\u001b[1;32m    503\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mx\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mndim\u001b[0m \u001b[0;34m>\u001b[0m \u001b[0;36m2\u001b[0m \u001b[0;32mor\u001b[0m \u001b[0my\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mndim\u001b[0m \u001b[0;34m>\u001b[0m \u001b[0;36m2\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mValueError\u001b[0m: x and y must have same first dimension, but have shapes (144,) and (100000,)"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAANQklEQVR4nO3cX4il9X3H8fenuxEak0aJk5DurmRb1pi90KITI6VpTUObXXuxBLxQQ6QSWKQx5FIpNLnwprkohKBmWWSR3GQvGkk2ZRMplMSCNd1Z8N8qynSlOl3BNYYUDFRWv704p51hnHWenXNmZp3v+wUD85znNzPf+TH73mfPznlSVUiStr7f2ewBJEkbw+BLUhMGX5KaMPiS1ITBl6QmDL4kNbFq8JMcSfJakmfPcz5JvptkPsnTSa6b/piSpEkNucJ/GNj3Huf3A3vGbweB700+liRp2lYNflU9BrzxHksOAN+vkSeAy5J8YloDSpKmY/sUPscO4JUlxwvjx15dvjDJQUb/CuDSSy+9/uqrr57Cl5ekPk6ePPl6Vc2s5WOnEfys8NiK92uoqsPAYYDZ2dmam5ubwpeXpD6S/OdaP3Yav6WzAOxacrwTODOFzytJmqJpBP8YcMf4t3VuBH5TVe96OkeStLlWfUonyQ+Am4ArkiwA3wI+AFBVh4DjwM3APPBb4M71GlaStHarBr+qblvlfAFfm9pEkqR14SttJakJgy9JTRh8SWrC4EtSEwZfkpow+JLUhMGXpCYMviQ1YfAlqQmDL0lNGHxJasLgS1ITBl+SmjD4ktSEwZekJgy+JDVh8CWpCYMvSU0YfElqwuBLUhMGX5KaMPiS1ITBl6QmDL4kNWHwJakJgy9JTRh8SWrC4EtSEwZfkpow+JLUhMGXpCYMviQ1YfAlqQmDL0lNGHxJamJQ8JPsS/JCkvkk965w/iNJfpLkqSSnktw5/VElSZNYNfhJtgEPAPuBvcBtSfYuW/Y14Lmquha4CfiHJJdMeVZJ0gSGXOHfAMxX1emqegs4ChxYtqaADycJ8CHgDeDcVCeVJE1kSPB3AK8sOV4YP7bU/cCngTPAM8A3quqd5Z8oycEkc0nmzp49u8aRJUlrMST4WeGxWnb8ReBJ4PeBPwLuT/J77/qgqsNVNVtVszMzMxc4qiRpEkOCvwDsWnK8k9GV/FJ3Ao/UyDzwEnD1dEaUJE3DkOCfAPYk2T3+j9hbgWPL1rwMfAEgyceBTwGnpzmoJGky21dbUFXnktwNPApsA45U1akkd43PHwLuAx5O8gyjp4DuqarX13FuSdIFWjX4AFV1HDi+7LFDS94/A/zldEeTJE2Tr7SVpCYMviQ1YfAlqQmDL0lNGHxJasLgS1ITBl+SmjD4ktSEwZekJgy+JDVh8CWpCYMvSU0YfElqwuBLUhMGX5KaMPiS1ITBl6QmDL4kNWHwJakJgy9JTRh8SWrC4EtSEwZfkpow+JLUhMGXpCYMviQ1YfAlqQmDL0lNGHxJasLgS1ITBl+SmjD4ktSEwZekJgy+JDUxKPhJ9iV5Icl8knvPs+amJE8mOZXkF9MdU5I0qe2rLUiyDXgA+AtgATiR5FhVPbdkzWXAg8C+qno5ycfWaV5J0hoNucK/AZivqtNV9RZwFDiwbM3twCNV9TJAVb023TElSZMaEvwdwCtLjhfGjy11FXB5kp8nOZnkjpU+UZKDSeaSzJ09e3ZtE0uS1mRI8LPCY7XseDtwPfBXwBeBv0ty1bs+qOpwVc1W1ezMzMwFDytJWrtVn8NndEW/a8nxTuDMCmter6o3gTeTPAZcC7w4lSklSRMbcoV/AtiTZHeSS4BbgWPL1vwY+FyS7Uk+CHwWeH66o0qSJrHqFX5VnUtyN/AosA04UlWnktw1Pn+oqp5P8jPgaeAd4KGqenY9B5ckXZhULX86fmPMzs7W3NzcpnxtSXq/SnKyqmbX8rG+0laSmjD4ktSEwZekJgy+JDVh8CWpCYMvSU0YfElqwuBLUhMGX5KaMPiS1ITBl6QmDL4kNWHwJakJgy9JTRh8SWrC4EtSEwZfkpow+JLUhMGXpCYMviQ1YfAlqQmDL0lNGHxJasLgS1ITBl+SmjD4ktSEwZekJgy+JDVh8CWpCYMvSU0YfElqwuBLUhMGX5KaMPiS1ITBl6QmBgU/yb4kLySZT3Lve6z7TJK3k9wyvRElSdOwavCTbAMeAPYDe4Hbkuw9z7pvA49Oe0hJ0uSGXOHfAMxX1emqegs4ChxYYd3XgR8Cr01xPknSlAwJ/g7glSXHC+PH/l+SHcCXgEPv9YmSHEwyl2Tu7NmzFzqrJGkCQ4KfFR6rZcffAe6pqrff6xNV1eGqmq2q2ZmZmYEjSpKmYfuANQvAriXHO4Ezy9bMAkeTAFwB3JzkXFX9aBpDSpImNyT4J4A9SXYD/wXcCty+dEFV7f6/95M8DPyTsZeki8uqwa+qc0nuZvTbN9uAI1V1Ksld4/Pv+by9JOniMOQKn6o6Dhxf9tiKoa+qv558LEnStPlKW0lqwuBLUhMGX5KaMPiS1ITBl6QmDL4kNWHwJakJgy9JTRh8SWrC4EtSEwZfkpow+JLUhMGXpCYMviQ1YfAlqQmDL0lNGHxJasLgS1ITBl+SmjD4ktSEwZekJgy+JDVh8CWpCYMvSU0YfElqwuBLUhMGX5KaMPiS1ITBl6QmDL4kNWHwJakJgy9JTRh8SWrC4EtSE4OCn2RfkheSzCe5d4XzX07y9Pjt8STXTn9USdIkVg1+km3AA8B+YC9wW5K9y5a9BPxZVV0D3AccnvagkqTJDLnCvwGYr6rTVfUWcBQ4sHRBVT1eVb8eHz4B7JzumJKkSQ0J/g7glSXHC+PHzuerwE9XOpHkYJK5JHNnz54dPqUkaWJDgp8VHqsVFyafZxT8e1Y6X1WHq2q2qmZnZmaGTylJmtj2AWsWgF1LjncCZ5YvSnIN8BCwv6p+NZ3xJEnTMuQK/wSwJ8nuJJcAtwLHli5IciXwCPCVqnpx+mNKkia16hV+VZ1LcjfwKLANOFJVp5LcNT5/CPgm8FHgwSQA56pqdv3GliRdqFSt+HT8upudna25ublN+dqS9H6V5ORaL6h9pa0kNWHwJakJgy9JTRh8SWrC4EtSEwZfkpow+JLUhMGXpCYMviQ1YfAlqQmDL0lNGHxJasLgS1ITBl+SmjD4ktSEwZekJgy+JDVh8CWpCYMvSU0YfElqwuBLUhMGX5KaMPiS1ITBl6QmDL4kNWHwJakJgy9JTRh8SWrC4EtSEwZfkpow+JLUhMGXpCYMviQ1YfAlqQmDL0lNDAp+kn1JXkgyn+TeFc4nyXfH559Oct30R5UkTWLV4CfZBjwA7Af2Arcl2bts2X5gz/jtIPC9Kc8pSZrQkCv8G4D5qjpdVW8BR4EDy9YcAL5fI08AlyX5xJRnlSRNYPuANTuAV5YcLwCfHbBmB/Dq0kVJDjL6FwDA/yR59oKm3bquAF7f7CEuEu7FIvdikXux6FNr/cAhwc8Kj9Ua1lBVh4HDAEnmqmp2wNff8tyLRe7FIvdikXuxKMncWj92yFM6C8CuJcc7gTNrWCNJ2kRDgn8C2JNkd5JLgFuBY8vWHAPuGP+2zo3Ab6rq1eWfSJK0eVZ9SqeqziW5G3gU2AYcqapTSe4anz8EHAduBuaB3wJ3Dvjah9c89dbjXixyLxa5F4vci0Vr3otUveupdknSFuQrbSWpCYMvSU2se/C9LcOiAXvx5fEePJ3k8STXbsacG2G1vViy7jNJ3k5yy0bOt5GG7EWSm5I8meRUkl9s9IwbZcCfkY8k+UmSp8Z7MeT/C993khxJ8tr5Xqu05m5W1bq9MfpP3v8A/gC4BHgK2Ltszc3ATxn9Lv+NwC/Xc6bNehu4F38MXD5+f3/nvViy7l8Y/VLALZs99yb+XFwGPAdcOT7+2GbPvYl78bfAt8fvzwBvAJds9uzrsBd/ClwHPHue82vq5npf4XtbhkWr7kVVPV5Vvx4fPsHo9Qxb0ZCfC4CvAz8EXtvI4TbYkL24HXikql4GqKqtuh9D9qKADycJ8CFGwT+3sWOuv6p6jNH3dj5r6uZ6B/98t1y40DVbwYV+n19l9Df4VrTqXiTZAXwJOLSBc22GIT8XVwGXJ/l5kpNJ7tiw6TbWkL24H/g0oxd2PgN8o6re2ZjxLipr6uaQWytMYmq3ZdgCBn+fST7PKPh/sq4TbZ4he/Ed4J6qent0MbdlDdmL7cD1wBeA3wX+LckTVfXieg+3wYbsxReBJ4E/B/4Q+Ock/1pV/73Os11s1tTN9Q6+t2VYNOj7THIN8BCwv6p+tUGzbbQhezELHB3H/grg5iTnqupHGzLhxhn6Z+T1qnoTeDPJY8C1wFYL/pC9uBP4+xo9kT2f5CXgauDfN2bEi8aaurneT+l4W4ZFq+5FkiuBR4CvbMGrt6VW3Yuq2l1Vn6yqTwL/CPzNFow9DPsz8mPgc0m2J/kgo7vVPr/Bc26EIXvxMqN/6ZDk44zuHHl6Q6e8OKypm+t6hV/rd1uG952Be/FN4KPAg+Mr23O1Be8QOHAvWhiyF1X1fJKfAU8D7wAPVdWWu7X4wJ+L+4CHkzzD6GmNe6pqy902OckPgJuAK5IsAN8CPgCTddNbK0hSE77SVpKaMPiS1ITBl6QmDL4kNWHwJakJgy9JTRh8SWrifwHXe3WluIZOawAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(n, sample_mean_cov_list, color='black')\n",
    "plt.xticks(range(1,simulations_number+1))\n",
    "plt.xlabel(\"Number of simulations\")\n",
    "plt.ylabel(\"Coefficient of variation of sample mean\")\n",
    "\n",
    "plt.xscale('log')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fourth_moment_list = []"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fourth_moment_calculation(realisations):\n",
    "    for i in range(simulations_number):\n",
    "        sample_fourth_moment=scipy.stats.moment(realisations[i], moment=4)\n",
    "\n",
    "        fourth_moment_list.append(sample_fourth_moment)\n",
    "\n",
    "        #print(second_moment_list)\n",
    "    #return fourth_moment_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fourth_moment_calculation(realisations_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sample_variance_cov_list = []\n",
    "sample_variance = np.var(r_array)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sample_variance_cov_calculation(fourth_moments):\n",
    "    for i in range(simulations_number):\n",
    "        sample_variance_cov = fourth_moments[i]/sample_variance\n",
    "\n",
    "        sample_variance_cov_list.append(sample_variance_cov)\n",
    "\n",
    "    #return sample_mean_cov_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sample_variance_cov_calculation(fourth_moment_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(n, sample_variance_cov_list, color='black')\n",
    "plt.xticks(range(1,simulations_number+1))\n",
    "plt.xlabel(\"Number of simulations\")\n",
    "plt.ylabel(\"Coefficient of variation of sample variance\")\n",
    "\n",
    "plt.xscale('log')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "masonry_unit_weight_realisations = r\n",
    "convergence_simulation_number = 100_000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%store masonry_unit_weight_realisations\n",
    "%store convergence_simulation_number"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.9.7"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "ad2bdc8ecc057115af97d19610ffacc2b4e99fae6737bb82f5d7fb13d2f2c186"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}