z0_from_TLS-CKEY_max.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "ef5ed08a",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "960cc6c4",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_csv(\"CKEY_TLS_10x10.xyzch\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "2e81bfe1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>X</th>\n",
       "      <th>Y</th>\n",
       "      <th>Z</th>\n",
       "      <th>Classification</th>\n",
       "      <th>HeightAboveGround</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>6305075</th>\n",
       "      <td>236361.987</td>\n",
       "      <td>172960.323</td>\n",
       "      <td>-1.263</td>\n",
       "      <td>2.0</td>\n",
       "      <td>0.000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2673899</th>\n",
       "      <td>236361.033</td>\n",
       "      <td>172961.005</td>\n",
       "      <td>1.593</td>\n",
       "      <td>1.0</td>\n",
       "      <td>2.940</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>16465826</th>\n",
       "      <td>236361.345</td>\n",
       "      <td>172960.145</td>\n",
       "      <td>2.830</td>\n",
       "      <td>1.0</td>\n",
       "      <td>3.724</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13420580</th>\n",
       "      <td>236360.613</td>\n",
       "      <td>172959.976</td>\n",
       "      <td>3.300</td>\n",
       "      <td>1.0</td>\n",
       "      <td>4.307</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>28616257</th>\n",
       "      <td>236357.635</td>\n",
       "      <td>172956.193</td>\n",
       "      <td>-1.398</td>\n",
       "      <td>2.0</td>\n",
       "      <td>0.000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>32872463</th>\n",
       "      <td>236355.457</td>\n",
       "      <td>172955.858</td>\n",
       "      <td>3.127</td>\n",
       "      <td>1.0</td>\n",
       "      <td>4.468</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2203974</th>\n",
       "      <td>236361.118</td>\n",
       "      <td>172960.993</td>\n",
       "      <td>1.505</td>\n",
       "      <td>1.0</td>\n",
       "      <td>2.847</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>34289781</th>\n",
       "      <td>236353.583</td>\n",
       "      <td>172957.128</td>\n",
       "      <td>1.922</td>\n",
       "      <td>1.0</td>\n",
       "      <td>3.281</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>37247834</th>\n",
       "      <td>236355.394</td>\n",
       "      <td>172957.353</td>\n",
       "      <td>2.838</td>\n",
       "      <td>1.0</td>\n",
       "      <td>4.005</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3854838</th>\n",
       "      <td>236360.784</td>\n",
       "      <td>172957.953</td>\n",
       "      <td>-0.823</td>\n",
       "      <td>2.0</td>\n",
       "      <td>0.000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                   X           Y      Z  Classification  HeightAboveGround\n",
       "6305075   236361.987  172960.323 -1.263             2.0              0.000\n",
       "2673899   236361.033  172961.005  1.593             1.0              2.940\n",
       "16465826  236361.345  172960.145  2.830             1.0              3.724\n",
       "13420580  236360.613  172959.976  3.300             1.0              4.307\n",
       "28616257  236357.635  172956.193 -1.398             2.0              0.000\n",
       "32872463  236355.457  172955.858  3.127             1.0              4.468\n",
       "2203974   236361.118  172960.993  1.505             1.0              2.847\n",
       "34289781  236353.583  172957.128  1.922             1.0              3.281\n",
       "37247834  236355.394  172957.353  2.838             1.0              4.005\n",
       "3854838   236360.784  172957.953 -0.823             2.0              0.000"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.sample(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "d9b8d455",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "38463608"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.count()[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "c4280ac5",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = df.sample(30000000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "87049c99",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "30000000"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.count()[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "232c7650",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "236362.41 236352.41 172964.91 172954.91 8.982 -1.54\n"
     ]
    }
   ],
   "source": [
    "x_max = df.X.max()\n",
    "x_min = df.X.min()\n",
    "y_max = df.Y.max()\n",
    "y_min = df.Y.min()\n",
    "z_max = df.Z.max()\n",
    "z_min = df.Z.min()\n",
    "print(x_max, x_min, y_max, y_min, z_max, z_min)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "2aa24c92",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.01 2.032 +/- 0.013 >>> took 1.0 minutes\n",
      "0.12 2.196 +/- 0.008 >>> took 11.0 minutes\n",
      "0.23 2.209 +/- 0.004 >>> took 19.0 minutes\n",
      "0.34 2.217 +/- 0.004 >>> took 24.0 minutes\n",
      "0.45 2.22 +/- 0.004 >>> took 34.0 minutes\n",
      "0.56 2.223 +/- 0.003 >>> took 40.0 minutes\n",
      "0.67 2.225 +/- 0.002 >>> took 50.0 minutes\n",
      "0.78 2.226 +/- 0.004 >>> took 58.0 minutes\n",
      "0.89 2.229 +/- 0.002 >>> took 65.0 minutes\n",
      "1.0 2.231 +/- 0.0 >>> took 7.0 minutes\n",
      "H* = 8.597059999999999   S* = 570.9757644074962   A* = 1100.0   z0 = 2.2312331387077764\n"
     ]
    }
   ],
   "source": [
    "# c = 0\n",
    "site_L = x_max - x_min\n",
    "site_W = y_max - y_min\n",
    "site_A = site_L * site_W\n",
    "voxel_dim = 1.0\n",
    "\n",
    "replicates = []\n",
    "for frac in np.around(np.linspace(0.01,1,10),2):\n",
    "    start = time.time()\n",
    "    z0_list = []\n",
    "    for rep in range(10):\n",
    "        window_max_h = []\n",
    "        voxel_count = []\n",
    "        voxels = 0\n",
    "        df_s = df.sample(frac=frac, replace=False)\n",
    "        pts = df_s.count()[0]\n",
    "        for i in np.arange(x_min, x_max, voxel_dim):\n",
    "            for j in np.arange(y_min, y_max, voxel_dim):\n",
    "                df_window = df_s[(df_s['X'].between(i, i+voxel_dim)) &\n",
    "                                 (df_s['Y'].between(j, j+voxel_dim))]\n",
    "                window_max_h.append(df_window.HeightAboveGround.max())\n",
    "                for k in np.arange(z_min, z_max, voxel_dim):\n",
    "                    df_voxel = df_s[(df_s['X'].between(i, i+voxel_dim)) &\n",
    "                                    (df_s['Y'].between(j, j+voxel_dim)) &\n",
    "                                    (df_s['Z'].between(k, k+voxel_dim))]\n",
    "                    voxel_count.append(df_voxel.count()[0])\n",
    "                    voxels += 1\n",
    "        mu = np.median([c for c in voxel_count if c != 0])\n",
    "        vi = [min(p/mu,1.0)*voxel_dim**3 for p in voxel_count]\n",
    "        S_star = np.sum(vi)\n",
    "        A_star = voxels * voxel_dim**3\n",
    "        H_star = np.mean(window_max_h)\n",
    "        z0_list.append(0.5*H_star*S_star/A_star)\n",
    "        if frac == 1.0: break\n",
    "    replicates.append([frac, z0_list, pts, pts/site_A, voxels])\n",
    "    end = time.time()\n",
    "    print(f\"{frac} {np.around(np.mean(z0_list),3)} +/- {np.around(np.std(z0_list),3)} >>> took {np.around((end - start)/60.0,0)} minutes\")\n",
    "\n",
    "# print results using densest (last) point cloud\n",
    "print(f\"H* = {H_star}   S* = {S_star}   A* = {A_star}   z0 = {0.5*H_star*S_star/A_star}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "46e554cc",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiwAAAGdCAYAAAAxCSikAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABAVklEQVR4nO3deVyVZf7/8TeLAhrgriQGqCVumfuCmpZh2aJjlpNp2qRmoblkk45lyzdlWiyzcUnHpUahJsUlt7RSxESbCGxMNFMpU0mtBNxQ4Pr94Y8zHlkPnCO38Ho+HudR576v6z6fy+1+c9/XdR83Y4wRAACAhbmXdQEAAABFIbAAAADLI7AAAADLI7AAAADLI7AAAADLI7AAAADLI7AAAADLI7AAAADL8yzrApwlJydHx44dk6+vr9zc3Mq6HAAAUAzGGGVkZOjGG2+Uu3vB11HKTWA5duyYGjRoUNZlAACAEjhy5IgCAwML3F9uAouvr6+kywP28/Mr42oAAEBxpKenq0GDBrbzeEHKTWDJvQ3k5+dHYAEA4DpT1HQOJt0CAADLI7AAAADLI7AAAADLI7AAAADLI7AAAADLI7AAAADLI7AAAADLI7AAAADLI7AAAADLI7AAAADLI7AAAADLI7AAAADLKzdffggAuL6cO3dO+/bty7P9/PnzSklJUXBwsHx8fPLsDw0NVZUqVa5FibAQAgsAoEzs27dPbdu2dbhfQkKC2rRp44KKYGUEFgBAmQgNDVVCQkKe7cnJyRo8eLCWLl2qpk2b5tvPKkp6lUjiSpGjCCwAgDJRpUqVQq+UNG3a1PJXUkp6lUgqmytFp06d0pgxY3T48GG77e7ZF3RzdalRo0bFPtbBgwd14A8px8Pbbvsdd9yh6dOnO6XeKxFYAJRYdna24uLidPz4cQUEBKhbt27y8PAo67IAp0tNTc33SsqFCxf0/vvv59n+888/a9q0aZoyZYpuuummfI954sQJbd261fbe09NTHTp0UOXKlZ1W99X+/ve/66OPPsqzvXU9d31w7w2SDuftVJC6Upv3zygxNcdu865du9SxY0f17du3lNXaI7AAKJGYmBg9++yzSklJsW0LDg7WjBkz1L9//7IrDHCBIUOG6PPPP3e437Rp0xxqv3jxYg0bNszhzymuSZMm6ejRo/leYRm6w/ErLN6BUsegvFdYnB1WJAILgBKIiYnRgAEDdN999yk6OlotWrTQnj17NH36dA0YMEDLly8ntECSlJGRoejoaGVnZxe7z88//yxJ+uSTT7Rr164i2+fk5Kh9+/bq0KFDiessylNPPaW9e/fq2LFjkiRvb2/16dNH99xzj8LCwlSpUiW79ufPn9ehQ4fUsGHDPHNYfvvtN3322Wdav3693fh69Oihnj17umwMklSrVi1FR0e79DNcxc0YY8q6CGdIT0+Xv7+/0tLS5OfnV9blAOVWdna2GjdurJYtW2rVqlVyd//f45xycnLUr18/7dmzRwcOHOD2ELRhwwb16dPH9t7Ts+ifk40xys7OloeHh9zc3Apsl5WVZfv/O++8s0RXQBx1+vRpxcTE6MUXX9SxY8fk4yk9/XAvvfXmm8Xrn5amPn3u0YULmdp3KkdhPXpp6tSp6ty5c7F+bcqj4p6/CSwAHLJ161b17NlT8fHx6tSpU5798fHx6tKli7Zs2aIePXpc+wILUdCKDun6evbH9fT8kpycHC1cuFBTp07VH3/8obFjx2ry5MmqVq1agX1yx1dQvcYYrV+/Xn/961+1d+9eDR48WNOnT1eDBg1cNo6CglPreu769skbSnTMq+d/VKtWTbt27dItt9xSouNdr4p7/q6YcQ5AiR0/flyS1KJFi3z3527PbWcl19uKjs8//1w//fSTWrdubXclKzk5WYMGDXL4eFFRUXbLhA8dOqTTp09r2LBhdsd3Jnd3d40YMUKPPPKI3nrrLb355pu2ADNq1Kh8J5gWtnooMTFREydO1JdffqkePXroww8/LPHvaWmEhISoYcOG8nLP1qjEC3n2nz17Vt9//72aN2+uqlWr5nuM2k0ryeuPr5WZmSnp8tWbxMTEChdYiosrLAAccj1cYVm/fr3ezOcSfXZ2ts6dO5dvn3Pnzik5OVlNmzbN96f6KlWq5LnFNXr0aD344IPOKTofhd0OcaYffvhBN9988zX5rGPHjmnq1KlauHChWjW7RdtXLdINBZzQr/bBhx9o1qz3lFWtoaa98bbuvffea/ZrdPHiRWVnZxf4TJWrffvtt2rbtq1DQTc9PV2+vr7XbExWwRUWAC7RrVs3BQcHa/r06fnOYYmMjFRISIi6detWZjVu3LhRiYmJuvfee4vd5/fff1dycrKCgoJUo0aNItt/9tlnWrdunUsDy6FDh7Ry5Upt2rRJmzZtUu7Plz6e0swXnlbXrmF27VMOp+iFF1/Ua//3fwoOCbZt/+WXo5o6daq+O3ZB5///tI8OHTooPDxc999/vxo3buyyMVwtNTVVhw4dkiR1vrmmbljWp4ge/zNU0tCRVdXzo1M6ePCgLl265NIlwFe6Fp/DD9uFI7AAcIiHh4dmzJihAQMGqF+/fpo8ebJtlVBkZKTWrl2r5cuXl/mE26CgIC1btqzY7YuaN3G1Ll26lKa8YgkMDFT79u31xx9/6I8//tB//vMfGWN0T/tGGmmWSnFL7do3k9TnyRukX16XfrHfHv54ZbX/Z7a+OXpJN998s9q3b6/27dsrNDT0mvxE//PPP2vKlClaunSpmjVrprVr16rPXT2kUwccOs7JU6fU/HiUJkyYoPfee09///vf9eCDD1a4qxIVEYEFgMP69++v5cuX69lnn7U7cYeEhFhiSbO3t7dSUlK0YcMG3XPPPcXqU9RTV6/0xRdfaP/+/brttttKUWXR8vupvn79+vrNraaG7sg70TMnJ0cXL2aqcmWvfOekZFU7Lx39QQcOHNCBAwc0e/ZsSZev5ISEhDh/AJKWLFmixx9/3PY+ICBAf/rTn/Sf//xH//nPf/Ltc+nSJZ06dUq1atXKs1xYkmoF3GT7M/jQQw9JkubPn68RI0a4ZAyFKWgCdHJyst1/82OlidzXA+awACgxqz7pNjU1VUOHDtWmTZs0dOhQvfPOO6pevXqpj5uenq7nnntO8+fPt034vBYrU4YOHVqsX9dTp05pzZo1euCBB1SrVq0i2x88eFDBwcGaPXt2gRNDS+upp57SvHnzJF2+5eHr61tkn4sXL+rkyZOqXbt2kbdijh49Kknq27evVq1aVep6HZU7V6Uk+BLHy1jWDMCprqeltNLlpa9LlizR+PHjVaVKFc2bN08PPPBAiY+3ceNGjRgxQqdPn9abb76pkSNHumxlTUmVZKKnFV1P4+DLD0uPSbcAnKqkS4LL6qTj5uamxx9/XOHh4Ro1apT69u2rRx55RLNmzSrW1Ydcf/zxhyZMmKAlS5YoPDxc8+fPV1BQkAsrx/WksFuJYWFh+W5HyRBYABcrLz+BhYaGKiEhIc/25ORkDR48WEuXLrV7xseV/cpS/fr1tWbNGs2fP1/jR4/ShcO7tGLFcrmpeJM0Rw9+VMnJ+/Te268rYtxzlpjcWdJ5E1b68wQ4isACOMmRI0fk4eGhG2+80W67Mx9Wtnv3bjVp0kTe3t6F9CqdkydPFviskpo1a+bZlvvE0mrVquW7/+TJk3bvjTGqVq1aoU86La0nn3xSy5Yt0zPPPCNJ+vLLL7Vr1y61rueumN6npPk9in2sZd0ldb9Bbd6YrEX/+kh33323bd9DDz2k1q1bO7n6ohX1Z2rw4MH5brfaLRaCFxzBHBbASXr16qX4+HhNmzZNY8aMsU2SLOwf5cKuTEj/+4c5NTVVY8aM0fLly7Vo0SK7VRfOdPLkSdWpU8clx75SrVq1dOLECZddrcg9bq1atXTq1Cnbdh9PqddtN+Xb5+LFizpx4oTq1KmTZ6Ln8ePH9f2vl2zPMKlevbrOnj2rwYMHa+HChS4ZQ2Gut/lEBSnphFWrBS+UDnNYgGvs3LlzOnfunMaPH6/x48erefPmatmyZYHtf//9d0nShx9+WOCDyowx+vjjj/N8jqvUqlVLCxYs0NixY22f4+/vb5sLkt+Xsx04cEARERGaPXt2nqelnj9/XitXrtSSJUts2wICArR48WKX3loZM2aMoqOj5e3trapVq+qNN97QqFGjCp0ka5vouW5dvidDY4wWLVqkCRMmyMvLS15eXiqrn/fKy7yJgm4zFid4oeIhsABOMnDgQPn4+Cg+Pl7nz5/X999/r++//15hYWHy8vLK0z4jI0PS5UmdV37rbK6DBw/qp59+sr1v1aqV6tatq65du7psDG5ubho+fLiGDx+uEydOKCoqSu+//75mzpyp1cujte6DmWp61cmiayNf9fwi+vLJ5YpbVVu2btVbL72kb1LS1alTJ40cOVIDBgwo1rLW0nJ3d9epU6d01113af78+QoODi71Md3c3PTEE0+od+/eevLJJ7V+/XrLrRK63pSX4IVrg8ACOMnYsWM1duxYSVJmZqYiIyM1ffp0nT19Uu/PfEnNmzWza3/+woX//RR51ZyUZVHLNHt2vHwaBendOfMVHh5+zcYhXX4AWXJyspKSknTkyBFJUnjrYDXd9qS0zb6tj6SmUp7tPSX1fFRq8767fvzxR+3evVutW7d2+cPWpMvf8dOjRw/17dvX6VdyAgMDtXbtWq1bt+6aPtIeqOiYwwK4SE5OjsaOHauvls8p8dfPJ7R7R23v+4uTKyvYiRMn1KpVK6Wmptq21a5dW4899pga3RSg6tmn8vS5dOmSfv/tN9WoWTPfp5Ju3/erFixZposXL9q2vfPOOxo3bpxLxlBS19OzP4DyhDksQBnau3evhg8frvj4eI0ZNVwZg4bK94bih5bkffv06quv6tPIERo1Zq9eeeUVlz2J9EpxcXF2YcXLy0sZGRm2R7jn5/Lj4C+qcuXKBd4icXd3l4eHh7KzsyXJNk/GCkuEAVwfCCyAE128eFF///vfNW3aNAUHBys2Nlbdu3d3+DhNb7xNH3Z7UG+//bZefvllxcTEaMGCBbrzzjtdUPX/PPjggw5PJM29MhEfH2/5KxMFra6RWEoLWB2BBXCSpKQkDRkyRMnJyfrrX/+qqVOnlup5KZUqVdLzzz+v/v37a8SIEerVq5cef/xxvffee9fkakt5VJxn4lwvzzABKhoCC+AkEyZMUGZmphISEtSqVSunHffmm2/Wl19+qYULF2rkyJEKCwvTE0884bTjVyQFLaOVWEoLWB2BBXCSLVu2SJJmzZpVrPZZWVk6ffq0qlWrlu/zTfLj5uamCxculLjGiq6wZbQSS2kBKyOwwNKuxyd67t27t1jtzp49q//+979q2bJlsW/xhIWFqWPHjqUpr8R4jDqAssSyZlhaeX509/W2jLY8/14AKDssa0a5cL1+Q3B5xGPUAZQlAgvK3O+//64FCxbo0qVLxe5z7NgxSdJXX32lw4cP59m/fv36PNs6deqkXr16lbzQCo7HqAMoSwQWlLl169Zp0qRJkqS6desWq09uuPn444/zfbrqlXJycnTy5Endeuut2r17d+mKBQCUCQILylyXLl1Ur149ZWRkaOLEiRo7dmyRISR3AmhREzq//fZbPfXUUzp58qQeeughZ5cOALhGCCwoc40aNdK+ffs0depUPf/88/rggw80d+7cQr+VuKjlqadPn9aLL76oOXPmqHnz5tq+fXuZ3bYo6eoaiRU2AJCLVUKwlMTERD311FPatWuXhg0bpjfeeEO1a9cudn9jjKKiovTss8/q7NmzeuWVVzRmzJgir9i4UklX10issAFQ/hX3/E1ggeXk5ORowYIFGj16tOrV9NN/t8Somr9/sfrOePttLVu2TIG39dScBUsUGBjo4mqLVtJnyUhcYQFQ/hFYKrjr8YFruf773//q6aef1vbt2zXlL/fptQbbHD5G92XSQ8+8pqeffloeHh4uqBIA4Aw8h6WCK86XvOWnLG9BZGRk6JVXXtHMmTNt35/TM6yjdOoHh46Tlp6uVqkfauzYsVq8eLHmzp1bZk+HBQA4B1dYyqnCJnoW9cC1a32F5dFHH1VUVJTdtsDAQFWuXLnAPjk5Obp06ZIqVaokd3f3fNscOnTI7v2QIUP04Ycflr5gAIDTcIWlgitqFU3Tpk0tM5kzKCjI9v9//vOfFRwcXGSf1NRULVmyRMOGDVO9evUKbJedna25c+fqzJkz2rVrlzPKBQCUAa6wVDBW/P6aS5cu6cCBA2rWrFmx+zg6jlOnTikzM1P169cvTakAACfjCguuG5UqVXIorJRErVq1XHp8AIBr5X/zHwAAwEIILAAAwPIcCiyRkZFq3769fH19VadOHfXr10/79+8vtE9MTIzuuusu1a5dW35+furcubM+++yzPO1WrFihZs2aycvLS82aNdPKlSsdGwkAACi3HAossbGxioiI0M6dO7V582ZlZWUpPDxcZ8+eLbDPtm3bdNddd2n9+vVKSEhQz549df/99ysxMdHWJj4+XgMHDtSQIUO0e/duDRkyRA8//DCrOophypQp2r59u8uOf/ToUY0aNUppaWku+wwAAIpSqlVCJ0+eVJ06dRQbG6vu3bsXu1/z5s01cOBATZ06VZI0cOBApaena8OGDbY2d999t6pXr67o6OhiHbOirhLy8vKSp6en1q9fr9tvv73I9o6srjl69Kh69OihH3/8UTt37rTUw9esuNoJAOC4a7JKKPen7ho1ahS7T05OjjIyMuz6xMfHa/z48XbtevfurZkzZ5amvArj3Llz6tGjh0aMGFFkaDl8+LAkaf369YV+S/DRo0f1/PPPO7XOkijpNx1b4SsGAADOU+LAYozRhAkT1LVrV7Vo0aLY/WbMmKGzZ8/q4Ycftm1LTU1V3bp17drVrVtXqampBR4nMzNTmZmZtvfp6ekOVF9+3HTTTfrxxx8lSQsWLNCCBQuK1e/FF18s9md4eXmpTp06JaqvtIr6ioHBgwfnu50rLwBQvpQ4sIwePVrfffedQ/MnoqOj9fLLL2v16tV5ToBubm52740xebZdKTIyUq+88opjRZdD+/bt08WLF3Xu3Dndf//92rVrl9avX1/gLbpz585p//79atKkSb5XII4dO6YePXrowoUL2rp1qxo2bCgPD49CH5PvSqGhoUpISMizvThf4ggAKD9KFFjGjBmjNWvWaNu2bQoMDCxWn48//lhPPPGEPvnkE/Xq1ctuX7169fJcTTlx4kSeqy5Xmjx5siZMmGB7n56ergYNGjgwivLBw8NDPj4+8vHx0eeff677779fDz74oDZs2KBu3brlae/j46MuXbrke6xjx47pnnvukZubm3bu3KlGjRq5uvwiFfYVA2FhYde4GgBAWXEosBhjNGbMGK1cuVJbt25VSEhIsfpFR0frL3/5i6Kjo3Xvvffm2d+5c2dt3rzZbh7Lpk2bCjyxSpdvU3h5eTlSfrlXpUoVffrpp7r//vvV//579OW/56llMW/X/XH6tCL+8heZixe1ZUusJcIKAAC5HAosERERioqK0urVq+Xr62u7KuLv72+7LD958mQdPXrU9q240dHReuyxx/Tuu++qU6dOtj4+Pj7y9/eXJI0dO1bdu3fX66+/rr59+2r16tX6/PPPXbpct7yqUqWKHnnkEf2xd6ta7oiQdhSvX3VJK++WXvzldjVs2NClNQIA4CiHljUXNKdk8eLFGjZsmCRp2LBhSklJ0datWyVJPXr0UGxsbJ4+Q4cO1ZIlS2zvly9frhdeeEGHDh1So0aNNG3aNPXv37/YA6moy5qvFhUVpSFDhmjEsEc159Xxci9kHtDVVq1erUHPvKzHRz6tf/zjH4XOIQIAwBmKe/7m25rLkdywMnToUP3zn/+Uu7vj37zwz3/+UyNGjFBERITee+89QgsAwKX4tuYKJjesPPbYYyUOK5I0fPhwSdKIESMkidACALAEAks5EB0drSFDhmjIkCGlCiu5hg8fLmOMRo4cKTc3N82aNYvQAgAoUwSW69xHH32kwYMHa/DgwVq4cKE8PDycctwRI0bIGKMnn3xSbm5uevfddwktAIAywxyW61y1atXUq1cvffzxx04LK1d69913NW7cOO3evVu33nqr048PAKjYinv+Lt29A5S58+fPq2fPni4JK5J055132j4HAICywi2h69zFixc1evRozZkzp1jtc3JylJmZKS8vr2LNdblw4YKkgpe0AwBwLRBYyoH27dura9euxWr766+/KioqSoMGDSr0qw+udMMNN3A7CABQpggs1zlHpyB9++23ioqK0rPPPsu3GQMArhvMYQEAAJZHYAEAAJZHYAEAAJZHYAEAAJZHYAEAAJZHYAEAAJZHYAEAAJZHYAEAAJZHYAEAAJZHYAEAAJZHYAEAAJZHYAEAAJZHYAEAAJbHtzWXU+fOndO+ffvybE9OTrb779VCQ0NVpUoVl9YGAICjCCzl1L59+9S2bdsC9w8ePDjf7QkJCWrTpo2rygIAoEQILOVUaGioEhIS8mw/f/68UlJSFBwcLB8fn3z7AQBgNW7GGFPWRThDenq6/P39lZaWJj8/v7IuBwAAFENxz99MugUAAJZHYAEAAJZHYAEAAJZHYAEAAJZHYAEAAJZHYAEAAJZHYAEAAJZHYAEAAJZHYAEAAJZHYAEAAJZHYAEAAJZHYAEAAJZHYAEAAJZHYAEAAJZHYAEAAJZHYAEAAJZHYAEAAJZHYAEAAJZHYAEAAJZHYAEAAJZHYAEAAJZHYAEAAJZHYAEAAJZHYAEAAJZHYAEAAJZHYAEAAJZHYAEAAJZHYAEAAJZHYAEAAJZHYAEAAJZHYAEAAJZHYAEAAJZHYAEAAJZHYAEAAJZHYAEAAJZHYAEAAJZHYAEAAJZHYAEAAJZHYAEAAJZHYAEAAJbnUGCJjIxU+/bt5evrqzp16qhfv37av39/oX2OHz+uQYMGqUmTJnJ3d9e4cePytFmyZInc3NzyvC5cuODQYAAAQPnkUGCJjY1VRESEdu7cqc2bNysrK0vh4eE6e/ZsgX0yMzNVu3ZtTZkyRa1atSqwnZ+fn44fP2738vb2dqQ8AABQTnk60njjxo127xcvXqw6deooISFB3bt3z7dPcHCw3n33XUnSokWLCjy2m5ub6tWr50g5AACggijVHJa0tDRJUo0aNUpdyJkzZxQUFKTAwEDdd999SkxMLLR9Zmam0tPT7V4AAKB8KnFgMcZowoQJ6tq1q1q0aFGqIkJDQ7VkyRKtWbNG0dHR8vb2VlhYmA4cOFBgn8jISPn7+9teDRo0KFUNAADAutyMMaYkHSMiIrRu3Tpt375dgYGBxerTo0cP3XbbbZo5c2ah7XJyctSmTRt1795ds2bNyrdNZmamMjMzbe/T09PVoEEDpaWlyc/Pr9jjAAAAZSc9PV3+/v5Fnr8dmsOSa8yYMVqzZo22bdtW7LDiCHd3d7Vv377QKyxeXl7y8vJy+mcDAADrceiWkDFGo0ePVkxMjL788kuFhIS4pChjjJKSkhQQEOCS4wMAgOuLQ1dYIiIiFBUVpdWrV8vX11epqamSJH9/f/n4+EiSJk+erKNHj+rDDz+09UtKSpJ0eWLtyZMnlZSUpMqVK6tZs2aSpFdeeUWdOnXSzTffrPT0dM2aNUtJSUmaPXu2M8YIAACucw4Flrlz50q6PBflSosXL9awYcMkXX5Q3M8//2y3v3Xr1rb/T0hIUFRUlIKCgpSSkiJJOn36tEaOHKnU1FT5+/urdevW2rZtmzp06ODgcAAAQHlU4km3VlPcSTsAAMA6inv+5ruEAACA5RFYAACA5RFYAACA5RFYAACA5RFYAACA5RFYAACA5RFYAACA5RFYAACA5RFYAACA5RFYAACA5RFYAACA5RFYAACA5RFYAACA5RFYAACA5RFYAACA5RFYAACA5RFYAACA5RFYAACA5RFYAACA5RFYAACA5RFYAACA5RFYAACA5RFYAACA5RFYAACA5RFYAACA5RFYAACA5RFYAACA5RFYAACA5RFYAACA5RFYAACA5RFYAACA5RFYAACA5RFYAACA5RFYAACA5RFYAACA5RFYAACA5RFYAACA5RFYAACA5RFYAACA5RFYAACA5RFYAACA5RFYAACA5RFYAACA5RFYAACA5RFYAACA5RFYAACA5RFYAACA5RFYAACA5RFYAACA5RFYAACA5RFYAACA5RFYAACA5RFYAACA5RFYAACA5RFYAACA5RFYAACA5RFYAACA5RFYAACA5RFYAACA5RFYAACA5RFYAACA5RFYAACA5RFYAACA5RFYAACA5TkUWCIjI9W+fXv5+vqqTp066tevn/bv319on+PHj2vQoEFq0qSJ3N3dNW7cuHzbrVixQs2aNZOXl5eaNWumlStXOlIaAAAoxxwKLLGxsYqIiNDOnTu1efNmZWVlKTw8XGfPni2wT2ZmpmrXrq0pU6aoVatW+baJj4/XwIEDNWTIEO3evVtDhgzRww8/rF27djk2GgAAUC65GWNMSTufPHlSderUUWxsrLp3715k+x49eui2227TzJkz7bYPHDhQ6enp2rBhg23b3XffrerVqys6OrpYtaSnp8vf319paWny8/NzaBwAAKBsFPf8Xao5LGlpaZKkGjVqlOYwio+PV3h4uN223r17a8eOHQX2yczMVHp6ut0LAACUTyUOLMYYTZgwQV27dlWLFi1KVURqaqrq1q1rt61u3bpKTU0tsE9kZKT8/f1trwYNGpSqBgAAYF0lDiyjR4/Wd999V+xbNkVxc3Oze2+MybPtSpMnT1ZaWprtdeTIEafUAQAArMezJJ3GjBmjNWvWaNu2bQoMDCx1EfXq1ctzNeXEiRN5rrpcycvLS15eXqX+bAAAYH0OXWExxmj06NGKiYnRl19+qZCQEKcU0blzZ23evNlu26ZNm9SlSxenHB8AAFzfHLrCEhERoaioKK1evVq+vr62qyL+/v7y8fGRdPlWzdGjR/Xhhx/a+iUlJUmSzpw5o5MnTyopKUmVK1dWs2bNJEljx45V9+7d9frrr6tv375avXq1Pv/8c23fvt0ZYwQAANc5h5Y1FzSnZPHixRo2bJgkadiwYUpJSdHWrVsL7RcUFKSUlBTb++XLl+uFF17QoUOH1KhRI02bNk39+/cvbmksawYA4DpU3PN3qZ7DYiUEFgAArj/X5DksAAAA1wKBBQAAWB6BBQAAWB6BBQAAWB6BBQAAWB6BBQAAWB6BBQAAWB6BBQAAWB6BBQAAWB6BBQAAWB6BBQAAWB6BBQAAWB6BBQAAWB6BBQAAWB6BBQAAWB6BBQAAWB6BBQAAWB6BBQAAWB6BBQAAWB6BBQAAWB6BBQAAWB6BBQAAWB6BBQAAWB6BBQAAWB6BBQAAWB6BBQAAWB6BBQAAWB6BBQAAWB6BBQAAWB6BBQAAWB6BBQAAWB6BBQAAWB6BBQAAWB6BBQAAWB6BBQAAWB6BBQAAWB6BBQAAWB6BBQAAWB6BBQAAWB6BBQAAWB6BBQAAWB6BBQAAWB6BBQAAWB6BBQAAWB6BBQAAWB6BBQAAWB6BBQAAWB6BBQAAWB6BBQAAWB6BBQAAWB6BBQAAWB6BBQAAWB6BBQAAWB6BBQAAWJ5nWRdwvcjOzlZcXJyOHz+ugIAAdevWTR4eHmVdFgAAFQJXWIohJiZGjRs3Vs+ePTVo0CD17NlTjRs3VkxMTFmXBgBAhUBgKUJMTIwGDBigli1bKj4+XhkZGYqPj1fLli01YMAAQgsAANeAmzHGlHURzpCeni5/f3+lpaXJz8/PKcfMzs5W48aN1bJlS61atUru7v/Ldzk5OerXr5/27NmjAwcOcHsIAIASKO75mysshYiLi1NKSor+9re/2YUVSXJ3d9fkyZN1+PBhxcXFlVGFAABUDASWQhw/flyS1KJFi3z3527PbQcAAFyDwFKIgIAASdKePXvy3Z+7PbcdAABwDQJLIbp166bg4GBNnz5dOTk5dvtycnIUGRmpkJAQdevWrYwqBACgYiCwFMLDw0MzZszQ2rVr1a9fP7tVQv369dPatWv11ltvMeEWAAAX48FxRejfv7+WL1+uZ599Vl26dLFtDwkJ0fLly9W/f/8yrA4AgIrBoSsskZGRat++vXx9fVWnTh3169dP+/fvL7JfbGys2rZtK29vbzVs2FDz5s2z279kyRK5ubnleV24cMGx0bhI//799eOPP2rLli2KiorSli1bdODAAcIKAADXiENXWGJjYxUREaH27dsrKytLU6ZMUXh4uPbu3auqVavm2+fw4cPq06ePRowYoaVLl+qrr77S008/rdq1a+vBBx+0tfPz88sTfry9vUswJNfw8PBQjx49yroMAAAqJIcCy8aNG+3eL168WHXq1FFCQoK6d++eb5958+bppptu0syZMyVJTZs21TfffKO33nrLLrC4ubmpXr16DpYPAAAqglJNuk1LS5Mk1ahRo8A28fHxCg8Pt9vWu3dvffPNN7p06ZJt25kzZxQUFKTAwEDdd999SkxMLPSzMzMzlZ6ebvcCAADlU4kDizFGEyZMUNeuXQt8sJokpaamqm7dunbb6tatq6ysLJ06dUqSFBoaqiVLlmjNmjWKjo6Wt7e3wsLCdODAgQKPGxkZKX9/f9urQYMGJR0KAACwuBIHltGjR+u7775TdHR0kW3d3Nzs3ud+fVHu9k6dOmnw4MFq1aqVunXrpn//+9+65ZZb9N577xV4zMmTJystLc32OnLkSEmHAgAALK5Ey5rHjBmjNWvWaNu2bQoMDCy0bb169ZSammq37cSJE/L09FTNmjXz7ePu7q727dsXeoXFy8tLXl5ejhcPAACuOw5dYTHGaPTo0YqJidGXX36pkJCQIvt07txZmzdvttu2adMmtWvXTpUqVSrwc5KSknjkPQAAkORgYImIiNDSpUsVFRUlX19fpaamKjU1VefPn7e1mTx5sh577DHb+1GjRumnn37ShAkTlJycrEWLFmnhwoWaOHGirc0rr7yizz77TIcOHVJSUpKeeOIJJSUladSoUU4YIgAAuN45dEto7ty5kpTneSSLFy/WsGHDJF3+5uKff/7Zti8kJETr16/X+PHjNXv2bN14442aNWuW3ZLm06dPa+TIkUpNTZW/v79at26tbdu2qUOHDiUcFgAAKE/cTO4M2Otcenq6/P39lZaWJj8/v7IuBwAAFENxz998+SEAALA8vvzwKufOndO+ffvy3Xf+/HmlpKQoODhYPj4+efaHhoaqSpUqri4RAIAKh8BylX379qlt27Yl6puQkKA2bdo4uSIAAEBguUpoaKgSEhLy3ZecnKzBgwdr6dKlatq0ab59AQCA8xFYrlKlSpUir5I0bdqUKykAAFxDTLoFAACWR2ABAACWR2ABAACWR2ABAACWV2EDy9y5c9W8eXOdPXvWJce/cOGC2rRpoxkzZrjk+AAAVCQVdpXQhg0btHfvXt1www3q27evPDw8iuzzxx9/SJImTpyo6tWrF9guJydHq1atkiR5e3vr2WefdUrNAABUVBU2sIwePVqffvqpJGn16tW6++67i+xz4cIF23/PnTtXYLuNGzfa/n/y5MmlrBQAAFTYwBIeHi5jjObMmaOIiAjde++9Gj16dKF9ch/bX9gj+BctWqSNGzfqjTfe0HPPPeeK0gEAqHAqbGDJ9fTTT+vgwYMaO3asgoODdd999xXYtqiHyn3++ed68skn9eSTT2rixImuKBcAgArJzRhjyroIZyju11PnJzs7Ww899JA+++wzxcXFlegptt9//726dOmizp07a+3atfL0rPBZEACAIhX3/F1hVwldycPDQ0uXLlXz5s1133336ciRIw71T01NVZ8+fRQUFKR///vfhBUAAJyMM+v/V6VKFX366afq3rm9Jg29W0uXLpWb3IrVd8rwJ9TE74IWrd7s8NUdAABQNALLFX7//XfV80jTsu5p0vwexe63sIOkDtL+n5KkkFtcVR4AABUWgeX/+/XXX9WnTx/V8g9UxqD35XvDDcXue/7CeQ0fPlw7F41T7I4uCgwMdGGlAABUPEy61eXlyj179tTPP/+sXbt26aabbnL4848fP66OHTuqRo0aiouLk6+vr8PHAACgomHSbTHl5OToscce0549e7R27doShRVJCggI0Lp163To0CENHDhQWVlZTq4UAICKq8IHlueff14xMTGKjo5W27ZtS3Wsli1bavny5dq0aZOeeeYZlZOLVwAAlLkKHVjmzZunt956SzNnztQDDzzglGOGh4dr3rx5mjt3rt5++22nHBMAgIquwk663bJliyIiIvTMM8/omWeeceqxhw8froMHD+q5555TaGio7r33XqceHwCAiqbCBpZZs2YpJydHjz76qBISEorV58KFC0pJSVFwcLC8vb0Lbdu/f39FR0drxYoVBBYAAEqpwgaWVatWSZI6duzo0s9htRAAAKVXYQNLfHy8Dh06pKZNm9ptP3/+vFJSUvLtc/jwYb3wwgt67bXXFBISkmd/cHCwfHx87LZdfXwAAOA4nsNylW+//bbEq4USEhJK9MWJAABUVMU9f1fYKywFCQ0NLXBOS+7Vl/yupOT2BQAAzscVFgAAUGZ40i0AACg3CCwAAMDyCCwAAMDyCCwAAMDyCCwAAMDyCCwAAMDyCCwAAMDyCCwAAMDyCCwAAMDyCCwAAMDyCCwAAMDyCCwAAMDyCCwAAMDyCCwAAMDyCCwAAMDyCCwAAMDyCCwAAMDyPMu6gOtFdna24uLidPz4cQUEBKhbt27y8PAo67IAAKgQuMJSDDExMWrcuLF69uypQYMGqWfPnmrcuLFiYmLKujQAACoErrAUISYmRgMGDNC9996r5557Tj4+Pjp//rw2bNigAQMGaPny5erfv39ZlwkAQLnmZowxZV2EM6Snp8vf319paWny8/NzyjGzs7PVuHFj1apVSydPntRPP/1k2xcUFKTatWvrt99+04EDB7g9BABACRT3/M0toULExcUpJSVF33zzjW699VbFx8crIyND8fHxuvXWW/XNN9/o8OHDiouLK+tSAQAo1wgshTh69Kgk6Z577tGqVavUqVMn3XDDDerUqZNWrVqle+65x64dAABwDQJLIU6ePClJ6t+/v9zd7X+p3N3d1a9fP7t2AADANQgshahdu7akyxNvc3Jy7Pbl5ORo1apVdu0AAIBrEFgKUb9+fUnSxo0b1a9fP7s5LP369dPGjRvt2gEAANdglVAhrlwldOrUKaWkpNj2hYSEqGbNmqwSAgCgFIp7/uY5LIXw8PDQjBkzbM9hmThxou05LBs3btS6deu0fPlywgoAAC5GYClC//79tXz5cj377LNau3atbXtISAgPjQMA4BrhllAx8V1CAAA4H7eEnMzDw0M9evQo6zIAAKiQWCUEAAAsj8ACAAAsj8ACAAAsz6HAEhkZqfbt28vX11d16tRRv379tH///iL7xcbGqm3btvL29lbDhg01b968PG1WrFihZs2aycvLS82aNdPKlSsdKQ0AAJRjDgWW2NhYRUREaOfOndq8ebOysrIUHh6us2fPFtjn8OHD6tOnj7p166bExET97W9/0zPPPKMVK1bY2sTHx2vgwIEaMmSIdu/erSFDhujhhx/Wrl27Sj4yAABQbpRqWfPJkydVp04dxcbGqnv37vm2ef7557VmzRolJyfbto0aNUq7d+9WfHy8JGngwIFKT0/Xhg0bbG3uvvtuVa9eXdHR0cWqxdXLmgEAgPMV9/xdqjksaWlpkqQaNWoU2CY+Pl7h4eF223r37q1vvvlGly5dKrTNjh07CjxuZmam0tPT7V4AAKB8KnFgMcZowoQJ6tq1q1q0aFFgu9TUVNWtW9duW926dZWVlaVTp04V2iY1NbXA40ZGRsrf39/2atCgQUmHAgAALK7EgWX06NH67rvvinXLxs3Nze597l2oK7fn1+bqbVeaPHmy0tLSbK8jR444Uj4AALiOlOhJt2PGjNGaNWu0bds2BQYGFtq2Xr16ea6UnDhxQp6enqpZs2ahba6+6nIlLy8veXl52d7nhiBuDQEAcP3IPW8XNaXWocBijNGYMWO0cuVKbd26VSEhIUX26dy5sz799FO7bZs2bVK7du1UqVIlW5vNmzdr/Pjxdm26dOlS7NoyMjIkiVtDAABchzIyMuTv71/gfodWCT399NOKiorS6tWr1aRJE9t2f39/+fj4SLp8q+bo0aP68MMPJV1e1tyiRQs9+eSTGjFihOLj4zVq1ChFR0frwQcflCTt2LFD3bt317Rp09S3b1+tXr1aL7zwgrZv366OHTsWq7acnBwdO3ZMvr6+hd5KKo309HQ1aNBAR44cuW5XIpWHMUjlYxzlYQwS47CS8jAGqXyMozyMQbo24zDGKCMjQzfeeKPc3QuZqWIcICnf1+LFi21thg4dam6//Xa7flu3bjWtW7c2lStXNsHBwWbu3Ll5jv3JJ5+YJk2amEqVKpnQ0FCzYsUKR0q7JtLS0owkk5aWVtallFh5GIMx5WMc5WEMxjAOKykPYzCmfIyjPIzBGGuNw+FbQkVZsmRJnm233367vv3220L7DRgwQAMGDHCkHAAAUEHwXUIAAMDyCCwO8PLy0ksvvWS3Oul6Ux7GIJWPcZSHMUiMw0rKwxik8jGO8jAGyVrjKNWj+QEAAK4FrrAAAADLI7AAAADLI7AAAADLI7AAAADLq9CBZc6cOQoJCZG3t7fatm2ruLi4QtvHxsaqbdu28vb2VsOGDTVv3jy7/d9//70efPBBBQcHy83NTTNnznRh9Zc5Mobjx49r0KBBatKkidzd3TVu3Lg8bRYsWKBu3bqpevXqql69unr16qWvv/7ahSO4zJFxxMTE6K677lLt2rXl5+enzp0767PPPsvTpl27dqpWrZqqVq2q2267Tf/6179cPQyHxrF9+3aFhYWpZs2a8vHxUWhoqN55550C23/00Udyc3NTv379XFD5/zj69yLXV199JU9PT912221225csWSI3N7c8rwsXLrig+v9xZBxbt27Nt8Z9+/bZtTt9+rQiIiIUEBAgb29vNW3aVOvXr7fEGCQpMzNTU6ZMUVBQkLy8vNSoUSMtWrTItr9Hjx75jvPee+912RgcHcewYcPyrbF58+Z27WbOnKkmTZrIx8dHDRo00Pjx4y31Z0qSli1bplatWqlKlSoKCAjQ448/rt9++822/9KlS3r11VfVqFEjeXt7q1WrVtq4caNLx3Clbdu26f7779eNN94oNzc3rVq1qsg+RZ0HXaqMH1xXZj766CNTqVIls2DBArN3714zduxYU7VqVfPTTz/l2/7QoUOmSpUqZuzYsWbv3r1mwYIFplKlSmb58uW2Nl9//bWZOHGiiY6ONvXq1TPvvPOOpcZw+PBh88wzz5gPPvjA3HbbbWbs2LF52gwaNMjMnj3bJCYmmuTkZPP4448bf39/88svv1hmHGPHjjWvv/66+frrr80PP/xgJk+ebCpVqmS+/fZbW5stW7aYmJgYs3fvXvPjjz+amTNnGg8PD7Nx40bLjOPbb781UVFRZs+ePebw4cPmX//6l6lSpYp5//3387RNSUkx9evXN926dTN9+/a1zBhynT592jRs2NCEh4ebVq1a2e1bvHix8fPzM8ePH7d7uZKj49iyZYuRZPbv329XY1ZWlq1NZmamadeunenTp4/Zvn27SUlJMXFxcSYpKckSYzDGmAceeMB07NjRbN682Rw+fNjs2rXLfPXVV7b9v/32m9349uzZYzw8POyeVl7W4zh9+rRdjUeOHDE1atQwL730kq3N0qVLjZeXl1m2bJk5fPiw+eyzz0xAQIAZN26cZcYRFxdn3N3dzbvvvmsOHTpk4uLiTPPmzU2/fv1sbf7617+aG2+80axbt84cPHjQzJkzx3h7e9v9W+ZK69evN1OmTDErVqwwkszKlSsLbV+c86ArVdjA0qFDBzNq1Ci7baGhoWbSpEn5tv/rX/9qQkND7bY9+eSTplOnTvm2DwoKcnlgcXQMV7r99tvzDSxXy8rKMr6+vuaDDz4oaZlFKs04cjVr1sy88sorhbZp3bq1eeGFF0pUY3E4Yxx/+tOfzODBg+22ZWVlmbCwMPPPf/7TDB061KWBpaRjGDhwoHnhhRfMSy+9lG9g8ff3d3KlhXN0HLmB5Y8//ijwmHPnzjUNGzY0Fy9edGapBXJ0DBs2bDD+/v7mt99+K/ZnvPPOO8bX19ecOXOmVLUWprR/L1auXGnc3NxMSkqKbVtERIS544477NpNmDDBdO3atfQFF8DRcbz55pumYcOGdttmzZplAgMDbe8DAgLMP/7xD7s2ffv2NY8++qiTqi6+4gQWR8+DzlYhbwldvHhRCQkJCg8Pt9seHh6uHTt25NsnPj4+T/vevXvrm2++0aVLl1xWa0FKMoaSOHfunC5duqQaNWo47ZhXcsY4cnJylJGRUWCNxhh98cUX2r9/v7p3717qmvPjjHEkJiZqx44duv322+22v/rqq6pdu7aeeOIJp9Wbn5KOYfHixTp48KBeeumlAtucOXNGQUFBCgwM1H333afExESn1X210vxetG7dWgEBAbrzzju1ZcsWu31r1qxR586dFRERobp166pFixaaPn26srOzLTGGNWvWqF27dnrjjTdUv3593XLLLZo4caLOnz9f4OcsXLhQf/7zn1W1alWn1p/LGX8vFi5cqF69eikoKMi2rWvXrkpISLDdrj506JDWr1/vsltbJRlHly5d9Msvv2j9+vUyxujXX3/V8uXL7WrMzMyUt7e3XT8fHx9t377d+YNwgrI+Dzr0XULlxalTp5Sdna26devaba9bt65SU1Pz7ZOamppv+6ysLJ06dUoBAQEuqzc/JRlDSUyaNEn169dXr169nHbMKzljHDNmzNDZs2f18MMP221PS0tT/fr1lZmZKQ8PD82ZM0d33XWX02q/UmnGERgYqJMnTyorK0svv/yyhg8fbtv31VdfaeHChUpKSnJF2XZKMoYDBw5o0qRJiouLk6dn/v+chIaGasmSJWrZsqXS09P17rvvKiwsTLt379bNN99siXEEBARo/vz5atu2rTIzM/Wvf/1Ld955p7Zu3WoLuYcOHdKXX36pRx99VOvXr9eBAwcUERGhrKwsTZ06tczHcOjQIW3fvl3e3t5auXKlTp06paefflq///673TyWXF9//bX27NmjhQsXOrX2K5X27/fx48e1YcMGRUVF2W3/85//rJMnT6pr164yxigrK0tPPfWUJk2a5NT6c5VkHF26dNGyZcs0cOBAXbhwQVlZWXrggQf03nvv2dr07t1bb7/9trp3765GjRrpiy++0OrVq10Sgp2hrM+DFTKw5HJzc7N7b4zJs62o9vltv5YcHYMj3njjDUVHR2vr1q15fgpwtpKOIzo6Wi+//LJWr16tOnXq2O3z9fVVUlKSzpw5oy+++EITJkxQw4YN1aNHD2eWbqck44iLi9OZM2e0c+dOTZo0SY0bN9YjjzyijIwMDR48WAsWLFCtWrVcVvPVijuG7OxsDRo0SK+88opuueWWAo/XqVMnderUyfY+LCxMbdq00XvvvadZs2Y5r/CrOPJ70aRJEzVp0sT2vnPnzjpy5IjeeustW2DJyclRnTp1NH/+fHl4eKht27Y6duyY3nzzTacHlpKMIScnR25ublq2bJn8/f0lSW+//bYGDBig2bNny8fHx679woUL1aJFC3Xo0MEltV+ppH+/lyxZomrVquWZaL5161ZNmzZNc+bMUceOHfXjjz9q7NixCggI0IsvvujM0u04Mo69e/fqmWee0dSpU9W7d28dP35czz33nEaNGmULie+++65GjBih0NBQubm5qVGjRnr88ce1ePFil42htMryPFghA0utWrXk4eGRJxmfOHEiT3rMVa9evXzbe3p6qmbNmi6rtSAlGYMj3nrrLU2fPl2ff/65br311lIfryClGcfHH3+sJ554Qp988km+V4Dc3d3VuHFjSdJtt92m5ORkRUZGuiSwlGYcISEhkqSWLVvq119/1csvv6xHHnlEBw8eVEpKiu6//35b25ycHEmSp6en9u/fr0aNGpXZGDIyMvTNN98oMTFRo0ePttVnjJGnp6c2bdqkO+64I08/d3d3tW/fXgcOHHBa7aUZR0E6deqkpUuX2t4HBASoUqVK8vDwsG1r2rSpUlNTdfHiRVWuXLn0xf9/JRlDQECA6tevbwsrufUZY/TLL7/YXc06d+6cPvroI7366qtOqzk/pfm9MMZo0aJFGjJkSJ5f2xdffFFDhgyxXY1s2bKlzp49q5EjR2rKlClyd3fubIeSjCMyMlJhYWF67rnnJEm33nqrqlatqm7duum1115TQECAateurVWrVunChQv67bffdOONN2rSpEm2fxOspqzPgxVyDkvlypXVtm1bbd682W775s2b1aVLl3z7dO7cOU/7TZs2qV27dqpUqZLLai1IScZQXG+++ab+7//+Txs3blS7du1KdayilHQc0dHRGjZsmKKioop939oYo8zMzFLVWxBn/X5cWWNoaKj++9//KikpyfZ64IEH1LNnTyUlJalBgwZlOgY/P7889Y0aNUpNmjRRUlKSOnbsWOAYk5KSXHb52Fm/F4mJiXY1hoWF6ccff7SFRkn64YcfFBAQ4NSwIpVsDGFhYTp27JjOnDljV5+7u7sCAwPt2v773/9WZmamBg8e7NS6r1aa34vY2Fj9+OOP+c7dOnfuXJ5Q4uHhIXN5IUnpC79KScZRUI2S8tTo7e2t+vXrKysrSytWrFDfvn2dWL3zlPl58JpM7bWg3CVqCxcuNHv37jXjxo0zVatWtc1EnzRpkhkyZIitfe5yrvHjx5u9e/eahQsX5lnOlZmZaRITE01iYqIJCAgwEydONImJiebAgQOWGIMxxlZf27ZtzaBBg0xiYqL5/vvvbftff/11U7lyZbN8+XK7pYUZGRkuGUNJxhEVFWU8PT3N7Nmz7Wo8ffq0rc306dPNpk2bzMGDB01ycrKZMWOG8fT0NAsWLLDMOP7xj3+YNWvWmB9++MH88MMPZtGiRcbPz89MmTKlwM9w9SqhkvyZulJ+q4Refvlls3HjRnPw4EGTmJhoHn/8cePp6Wl27dplmXG88847ZuXKleaHH34we/bsMZMmTTKSzIoVK2xtfv75Z3PDDTeY0aNHm/3795u1a9eaOnXqmNdee80SY8jIyDCBgYFmwIAB5vvvvzexsbHm5ptvNsOHD89z7K5du5qBAwe6pO7SjiPX4MGDTceOHfM95ksvvWR8fX1NdHS0OXTokNm0aZNp1KiRefjhhy0zjsWLFxtPT08zZ84cc/DgQbN9+3bTrl0706FDB1ubnTt3mhUrVpiDBw+abdu2mTvuuMOEhIQUulrNmTIyMmznBEnm7bffNomJibal2iU5D7pShQ0sxhgze/ZsExQUZCpXrmzatGljYmNjbfuGDh1qbr/9drv2W7duNa1btzaVK1c2wcHBZu7cuXb7Dx8+bCTleV19nLIcQ371BQUF2fYHBQXl2+bKZyCU9Thuv/32fGscOnSorc2UKVNM48aNjbe3t6levbrp3Lmz+eijj1w6BkfHMWvWLNO8eXNTpUoV4+fnZ1q3bm3mzJljsrOzCzy+qwOLMY7/mbpSfoFl3Lhx5qabbjKVK1c2tWvXNuHh4WbHjh0uqv5/HBnH66+/bho1amT789K1a1ezbt26PMfcsWOH6dixo/Hy8jINGzY006ZNs3tWS1mOwRhjkpOTTa9evYyPj48JDAw0EyZMMOfOnbNrs3//fiPJbNq0yWV1X83RcZw+fdr4+PiY+fPn53u8S5cumZdfftn2e9agQQPz9NNPu/xE7+g4Zs2aZZo1a2Z8fHxMQECAefTRR+2eabV161bTtGlT4+XlZWrWrGmGDBlijh496tIxXCl3OX9B/5aW5DzoSm7GuOD6GQAAgBNVyDksAADg+kJgAQAAlkdgAQAAlkdgAQAAlkdgAQAAlkdgAQAAlkdgAQAAlkdgAQAAlkdgAQAAlkdgAQAAlkdgAQAAlkdgAQAAlvf/AG/WbFDqQMTGAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# check plot, not for publication\n",
    "fig, ax = plt.subplots()\n",
    "ax.boxplot([y[1] for y in replicates],\n",
    "           notch=True,  # notch shape\n",
    "           vert=True,  # vertical box alignment\n",
    "           labels=[x[0] for x in replicates]);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "098fa598",
   "metadata": {},
   "outputs": [],
   "source": [
    "# export results for plotting in another notebook\n",
    "with open('CKEY_replicates_max', 'w') as outfile:\n",
    "    outfile.write('\\n'.join(str(i) for i in replicates))\n",
    "    \n",
    "# once this was complete, commas were manually inserted in the last row (frac = 1.0)\n",
    "# so that fields r1 through r9, representing replicates 2 through 10, would match\n",
    "# the overall schema."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python3 (mangrove)",
   "language": "python",
   "name": "mangrove"
  },
  "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.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}