Implement mixed precision: f32 stencil, f64 pedestal

- Stencil arithmetic and shared memory use float (COMPUTE_TYPE alias). - Pedestal accumulation stays double to preserve variance accuracy. Notes: - On RTX 4090, FP32 throughput is ~64× higher than FP64, so moving stencil math to float improves performance. - Using float also avoids shared memory bank conflicts: stride-18 maps to distinct banks for 32-bit values, but caused conflicts with 64-bit.
2026-11-07 21:53:53 +01:00 · 2026-04-27 14:56:40 +02:00
1 parent 40f08fad92
commit ac96d1f688
@@ -110,7 +110,7 @@
            "metadata": {},
            "outputs": [],
            "source": [
-                "N_STREAMS  = 1 \n",
+                "N_STREAMS  = 10\n",
                "cf_cuda = ClusterFinderCUDA(image_size, cluster_size, capacity=capacity, n_streams=N_STREAMS)"
            ]
        },
@@ -124,7 +124,7 @@
                    "name": "stdout",
                    "output_type": "stream",
                    "text": [
-                        "Pedestal (1000 frames): 0.498s\n"
+                        "Pedestal (1000 frames): 0.480s\n"
                    ]
                }
            ],
@@ -155,7 +155,7 @@
                    "name": "stdout",
                    "output_type": "stream",
                    "text": [
-                        "Reading 40000 frames:        2.002s  (19985 FPS, 6098.804 GB/s)\n"
+                        "Reading 40000 frames:        1.922s  (20815 FPS, 6352.256 GB/s)\n"
                    ]
                }
            ],
@@ -186,7 +186,7 @@
                    "name": "stdout",
                    "output_type": "stream",
                    "text": [
-                        "CPU clustering:          73.978s (541 FPS, 55449596 clusters, 1386.24/frame)\n"
+                        "CPU clustering:          73.207s (546 FPS, 55449596 clusters, 1386.24/frame)\n"
                    ]
                }
            ],
@@ -230,32 +230,23 @@
            "execution_count": 10,
            "id": "4b8df93b-9a1b-41a5-9fed-9cda295fb523",
            "metadata": {},
-            "outputs": [
-                {
-                    "name": "stdout",
-                    "output_type": "stream",
-                    "text": [
-                        "CUDA clustering:          3.880s  (10310 FPS, 55810704 clusters, 1395.27/frame)\n",
-                        "Speedup (CPU / CUDA):     19.07×\n"
-                    ]
-                }
-            ],
+            "outputs": [],
            "source": [
-                "# Warmup: first kernel launch pays CUDA context + pedestal H2D upload cost\n",
-                "cf_cuda.find_clusters(data[0])\n",
-                "_ = cf_cuda.steal_clusters(realloc_same_capacity=False)\n",
+                "# # Warmup: first kernel launch pays CUDA context + pedestal H2D upload cost\n",
+                "# cf_cuda.find_clusters(data[0])\n",
+                "# _ = cf_cuda.steal_clusters(realloc_same_capacity=False)\n",
                "\n",
-                "t0 = time.perf_counter()\n",
-                "for frame in data:\n",
-                "    cf_cuda.find_clusters(frame)\n",
-                "t_cuda = time.perf_counter() - t0\n",
-                "clusters_cuda = cf_cuda.steal_clusters(realloc_same_capacity=False)\n",
-                "n_clusters_cuda = clusters_cuda.size\n",
-                "print(f'CUDA clustering:          {t_cuda:.3f}s  ({N/t_cuda:.0f} FPS, '\n",
-                "      f'{n_clusters_cuda} clusters, {n_clusters_cuda/N:.2f}/frame)')\n",
-                "print(f'Speedup (CPU / CUDA):     {t_cpu / t_cuda:.2f}×')\n",
+                "# t0 = time.perf_counter()\n",
+                "# for frame in data:\n",
+                "#     cf_cuda.find_clusters(frame)\n",
+                "# t_cuda = time.perf_counter() - t0\n",
+                "# clusters_cuda = cf_cuda.steal_clusters(realloc_same_capacity=False)\n",
+                "# n_clusters_cuda = clusters_cuda.size\n",
+                "# print(f'CUDA clustering:          {t_cuda:.3f}s  ({N/t_cuda:.0f} FPS, '\n",
+                "#       f'{n_clusters_cuda} clusters, {n_clusters_cuda/N:.2f}/frame)')\n",
+                "# print(f'Speedup (CPU / CUDA):     {t_cpu / t_cuda:.2f}×')\n",
                "\n",
-                "hist_cuda = make_hist(clusters_cuda)"
+                "# hist_cuda = make_hist(clusters_cuda)"
            ]
        },
        {
@@ -263,27 +254,36 @@
            "execution_count": 11,
            "id": "1fda99ee-db54-40c1-8e5b-51168a33fb96",
            "metadata": {},
-            "outputs": [],
+            "outputs": [
+                {
+                    "name": "stdout",
+                    "output_type": "stream",
+                    "text": [
+                        "CUDA clustering:          2.385s  (16774 FPS, 55081141 clusters, 1377.03/frame)\n",
+                        "Speedup (CPU / CUDA):     30.70×\n"
+                    ]
+                }
+            ],
            "source": [
-                "# BATCH_SIZE = 500\n",
+                "BATCH_SIZE = 2000\n",
                "\n",
-                "# # warmup\n",
-                "# _ = cf_cuda.find_clusters_batched(data[:1], first_frame=0)\n",
+                "# warmup\n",
+                "_ = cf_cuda.find_clusters_batched(data[0:BATCH_SIZE], first_frame=0)\n",
                "\n",
-                "# t0 = time.perf_counter()\n",
-                "# clusters_cuda_per_frame = []\n",
-                "# for start in range(0, N, BATCH_SIZE):\n",
-                "#     stop = min(start + BATCH_SIZE, N)\n",
-                "#     clusters_cuda_per_frame.extend(\n",
-                "#         cf_cuda.find_clusters_batched(data[start:stop], first_frame=start)\n",
-                "#     )\n",
-                "# t_cuda = time.perf_counter() - t0\n",
+                "t0 = time.perf_counter()\n",
+                "clusters_cuda_per_frame = []\n",
+                "for start in range(0, N, BATCH_SIZE):\n",
+                "    stop = min(start + BATCH_SIZE, N)\n",
+                "    clusters_cuda_per_frame.extend(\n",
+                "        cf_cuda.find_clusters_batched(data[start:stop], first_frame=start)\n",
+                "    )\n",
+                "t_cuda = time.perf_counter() - t0\n",
                "\n",
-                "# n_clusters_cuda = sum(cv.size for cv in clusters_cuda_per_frame)\n",
+                "n_clusters_cuda = sum(cv.size for cv in clusters_cuda_per_frame)\n",
                "\n",
-                "# print(f'CUDA clustering:          {t_cuda:.3f}s  ({N/t_cuda:.0f} FPS, '\n",
-                "#       f'{n_clusters_cuda} clusters, {n_clusters_cuda/N:.2f}/frame)')\n",
-                "# print(f'Speedup (CPU / CUDA):     {t_cpu / t_cuda:.2f}×')"
+                "print(f'CUDA clustering:          {t_cuda:.3f}s  ({N/t_cuda:.0f} FPS, '\n",
+                "      f'{n_clusters_cuda} clusters, {n_clusters_cuda/N:.2f}/frame)')\n",
+                "print(f'Speedup (CPU / CUDA):     {t_cpu / t_cuda:.2f}×')"
            ]
        },
        {
@@ -293,14 +293,14 @@
            "metadata": {},
            "outputs": [],
            "source": [
-                "# def make_hist_from_batch(result_list):\n",
-                "#     h = bh.Histogram(bh.axis.Regular(100, -2, 4000))\n",
-                "#     energies = [np.asarray(cv.sum()).ravel() for cv in result_list if cv.size > 0]\n",
-                "#     if energies:\n",
-                "#         h.fill(np.concatenate(energies))\n",
-                "#     return h\n",
+                "def make_hist_from_batch(result_list):\n",
+                "    h = bh.Histogram(bh.axis.Regular(100, -2, 4000))\n",
+                "    energies = [np.asarray(cv.sum()).ravel() for cv in result_list if cv.size > 0]\n",
+                "    if energies:\n",
+                "        h.fill(np.concatenate(energies))\n",
+                "    return h\n",
                "\n",
-                "# hist_cuda = make_hist_from_batch(clusters_cuda_per_frame)"
+                "hist_cuda = make_hist_from_batch(clusters_cuda_per_frame)"
            ]
        },
        {
@@ -323,7 +323,7 @@
                    "name": "stdout",
                    "output_type": "stream",
                    "text": [
-                        "Cluster count diff:       361108 (0.65%)\n"
+                        "Cluster count diff:       368455 (0.66%)\n"
                    ]
                }
            ],
@@ -350,7 +350,7 @@
                {
                    "data": {
                        "image/png":
-                            "iVBORw0KGgoAAAANSUhEUgAAA3kAAAJOCAYAAAAK+M50AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAlypJREFUeJzt3Qd4FNXawPF3Nx1SSGiBEJpUqQIqKKKogGDBjvcqomIHFVFR7Ni4lqvop1gRxIZXEUVABa4ConKVJr1EaQIhtDTSs/M978GN2TSSkM22/+95lsnOzM6ePXt2mXffc87YLMuyBAAAAADgF+yeLgAAAAAAoOYQ5AEAAACAHyHIAwAAAAA/QpAHAAAAAH6EIA8AAAAA/AhBHgAAAAD4EYI8AAAAAPAjBHkAAAAA4EcI8gAAAADAjxDkAfAra9askeuvv15atWol4eHhEhkZKT169JDnnntODh06VLTfWWedZW7uMnnyZJk2bZrbjg/f4Y9t4YcffpArr7xSEhISJDQ0VGJiYuS0006T119/XY4cOVK0X8uWLcVmsxXd9PN46qmnyvTp012Op/tdcMEFZT7X8uXLzWO9pQ4r+x1T1df0+OOPu9RVnTp1pFmzZjJo0CD5v//7P8nIyCi3TPn5+RIfH28e99lnn9XwKwbgiwjyAPiNt99+W3r27Cm//vqr3HffffLNN9/IrFmz5IorrpA33nhDRo4cWWtl8ccTe1SPv7WFxx57TPr16ye7d++WJ598UhYsWCAzZsyQc845xwQqDz/8sMv+p59+uvz888/mpvWggciIESNMQOhrauM7Ro+pdaXLF154QZo3by7jxo2TTp06yW+//VbmY+bMmSP79u0zf0+ZMuW4ywDA9wV7ugAAUBP0pOi2226TAQMGyBdffCFhYWFF23TdPffcY06afJllWZKTkyMRERHiC7Kzs32mrN5CMzIaBAUHe+d/z59++qk88cQTJpjRgEfL6jR48GATjOhnsbh69epJ7969i+6fe+650qJFC3nxxRfNZ9ZX1NZ3jAaRDRo0KLp/1VVXyejRo+XMM8+Uiy66SLZs2eLy3M7ATjOqus/8+fPlzz//NFlAAIGLTB4Av/DMM8+YE8633nqr1AmQ0hMgPUEqz6JFi8zjdVnc9u3bS3Wr+uOPP8yJV9OmTc1zNW7c2GQxVq9eXdRNa/369bJ48eKirle6zik9PV3uvfde091Ly6Vd3saMGePSzU3p4/TkTjMEHTt2NM/13nvvVVgPn3zyifTp00fq1q1rupFpV69Vq1a57HPdddeZbUlJSTJkyBDzd2JiojlJzc3Nddk3Ly9PnnrqKenQoYN5/oYNG5quavv373fZz9k17fPPP5eTTjrJdGObMGGC2aZ1MXDgQNP9TB8/atQomTt3rkt9a0ZIA5tdu3aVek033HCD1K9f3wS45TnWe1K8jJp56dq1qylj69at5ZVXXil1vMq+Rw6Hw3Sl6969uwlonQHN7Nmzj9kWnG3u/fffN3Wvz6Fl1/fF2XWvJGcmTNtlydel2Rytey2Hthe973yM3tc2ccopp5iugtWlAV5sbKyps7LKFxUVZd7rimgdtW/fXnbs2CE1RTNcWp6yslhff/212eZ8T7Tt3nzzzabNO9u0ZhsXLlzo1u+Y49GtWzd56KGHZOfOneYzXtyePXtMcHnhhRea7KK2SX/KHAOoHu/8qRAAqqCwsFC+++478wu4nri5mwZG+pw6Bke7Uh04cEB++uknSU1NNds1iLj88svNOCXtqqecJ4VZWVnm13b9pf3BBx80wYYGAY8++qisXbvWnGgWP3nWjIGOf9LtOuamUaNGFZ6Ealc5DcJ0qQHa888/L2eccYb88ssvcuKJJ7pkjPSEVDMyGmAsWbLEBFpaZn0upSeLQ4cONc+vGRodc6Un5tpdT8czarBQPFO3cuVK2bhxo3luDY40qNi7d695vfq3ds/T8n/88ccmeC3ulltukaefflrefPNNE1Q66Rgn7Qqo+2tQVt33xEmDPg3WNIjS+vzwww/lrrvuMnWlQV1V3yMNmD/44ANTjxoA6Ym+1oMzCKuoLTiNHz/eBOYazNvt9grf44qCHD2OBgL6XBpgX3rppWbdf//736IA5f777zcB4bZt24reOy2rvl/ahbKi4EDfy3Xr1smwYcNMwF5d2va0HWlwVZNBkAa4U6dOLdVlUl+T1qm2ETV8+HDzHml7a9eunWkjev/gwYNe8x1TFv286udQP6vXXnuty+vT8umPIc4s6bvvvmvaQlmBOIAAYQGAj0tOTrb06+yqq66q9GPOPPNMc3P6/vvvzTF0Wdy2bdvM+qlTp5r7Bw4cMPcnTZpU4fE7derkcnyniRMnWna73fr1119d1n/22WfmuPPmzStap/djYmKsQ4cOHfP17Ny50woODrbuuOMOl/UZGRlWfHy8deWVVxatGzFihDn2f/7zH5d9hwwZYrVv377o/scff2z2mzlzpst+WnZdP3ny5KJ1LVq0sIKCgqzNmze77HvfffdZNpvNWr9+vcv6QYMGlapvLVejRo2s3NzconXPPvusqS99H8pT2fdEy6hlWb16tcv6AQMGWNHR0daRI0eq9B4tWbLE3H/ooYeq1Racba5fv36ltj322GNmW0naDnV98frQ1xUREWH9+eefRev0Nep+TZo0KXpd6osvvjDrZ8+eXbRu+/bt5r274YYbKnwdy5YtM4994IEHrMrSsmm7ys/PNzctt7P9adsovt/5559f5jGc7c35GSzPK6+8YvYr3gb1sxMWFmbdc889ResiIyOtMWPGWO7+jqnqa3K+5/v37y/zMdnZ2Wb74MGDi9Y5HA6rTZs2VkJCglVQUOBynP/+979VeIUA/E1Ad9fUX8O0e4N279Ffu/QX86rS8zAdGK2/Buqvs/oLn/5iCsA/xcXFyQknnGAyZDqmSLtCasarsrQLXefOnU33voKCgqKbdqssq7vo2WefbbrHHcu3335rjqO/8Bc/rma/NCtV8rj6XPr9V5xmrIp3odOyatc63a/4MbXsmgUreUx9vH4XFqfdFPX1Fs8iqn/84x+lXoNm1FJSUsy4L6X1qtm/888/36W76/G8Jzp5hWZ9ivvnP/9pumdqNqcq75F2A1Ta/fR4XHbZZXK8tKza3dNJu2cqzbgWz7o51xd/nzXzo6/PXRN2zJs3T0JCQsxNM4b/+c9/5I477nDJ2NaEq6++2vw/XDwbqVlj7YKs2W0n7bKq++jzL1u2zGQWfcHR332k1OdLu/dqFjYoKMis09eq7VSzeQACV0AHeTq2Qv+zf/XVV6t9DD0peeedd0ygt2nTJvnqq6/MfyAAao9OUqAnstoFzd305Em7v+kJv3YN1KnTtdvZnXfeWeEU5046A55Owe486XXedCyTnsRpN8PimjRpUqlyOWfWO/nkk0sdW8fwlDyu1lfJ7o96glx83JseU7uyaRfEksdMTk6uVFm1C5yOjyuprHXa3U67lr722mtFwZZ2JSzZtfN43hMNTktyrnN216vse6Rju/TEuqxjVkVl3+OKaKBbnL5nFa2vaHxjebQbrKrq56xv375mNkrt3rthwwbTpnRMn7MsSsdjapfDsmgAqvQ9qIi+Vu3SqJdncB5Lgzn9P1mDeyf9PGhQpP93azdZfZz+OKJtuia/Y2riNRXnDMz1h2knZ2B+ySWXmHrVm3bX1TqfOXNmqe7KAAJHQI/J05nA9FYeHaOhY0t0zIZ+Ueovu88++2zRtbV07In+yqxjFHQQOQDP0BNtnWRDMyvVnVXOGfCUnHikZCDjzHw4T650pjvNTOgYL/3O0HFVFdGTRR0LVd6v7MVn1VOVHVPjfJxeI0vLVxP0mDrhSXkzBmrQc6yy6uOdAWhx5Z1Qa2Cm09FrVk1/gNPMoM5ceCyVfU/Kel7nOi1rVd4jDST1JF4ffzyBWln1Vrw9Fh/DV1Z7rC36Grt06WJmb9Rxi5Udl6dBR69evSrcR4N+vSRDWZzry/phoCTNYmkmWC/roEGpBpclL9Wg79+kSZPMTScy0QlZHnjgAZNFLq+tV+c7pqZek5Nz4hjnOUhaWpoJ5Jw/7pTlo48+kttvv73SzwHAfwR0Jq8y/1n8+OOPZtC//qqrJx7nnXeebN261WzXrJ3OzKa/NmsXFO1OdOONN7pcDBVA7dAJJjTLctNNN5kT+5K0S5Z+Zsvj7A6on/WyTqzKo0GI/hikJ7/O7n5KT8z1EgIl6aQXv//+uwko9MS35K2ibokV0SyWZg702GUd91gn2WXRsmp2SwOZso5XmR+3tKuo/hCmGZzi9Hu1LJqR0JNznQxGJzjRE9SqTh5R3nuidAKVktca0xNhDVg1A1iV98j5I+GxrvdWXluoSHntsaI2XBseeeQROXz4sAnGy+o+mJmZaYLAqtIJQ8pqJ0oDdudF1I9FZ/bUbqs6AYveNFguq2uwk7Y1zRTrDwkl28rxfsfU1GtS2mZ1KIi2C70IvbPdarvSCZO+//77UjcNZumyCQSugM7kVUT/g9e+/PqLnbNrhM68pr/y6X8c+mWrU3Zr9wn91dDZPeTuu+82M6npLFwAao92u9KTbQ0KdAY8vZ6VdtHSEy8do6XTnms2vuQ4NCftcqcnZRMnTjRj4DQzpF0A9ZIAxelJt54U6o8+bdu2NV3O9POu6zUb4KQBhgYy2jVMfwzSk01dpzM76q/vejFp/b7QcWw6fkwzCnpyrMFNZU/8itOTP53dUWfU0+8m/UFKX4dm0XRmTZ3d0nlJg8rSSxJoTwadlVC7pmu3N+1ept+LehKpM29qUFYRfb16oqkBkZZPMxd6cqrd25XOJlkyY6Jj3HQWSC2zzl55LJV9T5R+n2uXPs3yaWZKZ8bUrI/20nBmpir7HmnXUp2pUcd2aT1rcKgBnbY3PZaOO6uoLVRE61y7ETpn7dQAXrselnWJieOl/4/pmEbtwniscXlaxxLine truncated
+                            "iVBORw0KGgoAAAANSUhEUgAAA3kAAAJOCAYAAAAK+M50AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAnI9JREFUeJzt3Qd4VFX6x/FfegHSIAESqiKIFLGuHSsqNuzu2ntXLOvaFXV114r+FXvDhhUrNlYFXXVVEEFUBAXpJqGkQHrm/7wnTswkkxBCJtO+n+e5XObeOzNnztyZ3HfeU2I8Ho9HAAAAAICIEBvsAgAAAAAA2g9BHgAAAABEEII8AAAAAIggBHkAAAAAEEEI8gAAAAAgghDkAQAAAEAEIcgDAAAAgAhCkAcAAAAAEYQgDwAAAAAiCEEegIgye/ZsnXrqqerfv7+Sk5PVuXNnbbvttrr99tu1evXq+uP23HNPtwTKhAkT9NRTTwXs8RE+IvFc+PTTT3XMMccoLy9PiYmJSk9P1y677KIHH3xQ69atqz+uX79+iomJqV/s8/iXv/xFEydO9Hk8O+7ggw/2+1zffPONu2+o1GFrv2M29jXdeOONPnWVmpqqXr16af/999f//d//qaSkpNkyVVVVqUePHu5+r7zySju/YgDhiCAPQMR49NFHtd122+nrr7/W3//+d7333nuaPHmyjj76aD300EM6/fTTO6wskXhhj7aJtHPhhhtu0B577KFly5bp5ptv1ocffqhJkyZpn332cYHKtdde63P8rrvuqi+++MItVg8WiJx88skuIAw3HfEdY49pdWXrO++8U3369NEVV1yhIUOG6LvvvvN7n7ffflu///67+//jjz++yWUAEP7ig10AAGgPdlF07rnnar/99tPrr7+upKSk+n227bLLLnMXTeHM4/GovLxcKSkpCgdlZWVhU9ZQYRkZC4Li40Pzz/PLL7+sm266yQUzFvBYWb0OPPBAF4zYZ7GhjIwM7bTTTvW39913X/Xt21d33323+8yGi476jrEgslu3bvW3jzvuOF1wwQUaOXKkDj30UP38888+z+0N7Cyjasd88MEHWrp0qcsCAoheZPIARIRbb73VXXA+8sgjTS6AjF0A2QVScz755BN3f1s3tGjRoibNqn799Vd34ZWbm+ueq3v37i6LMWvWrPpmWnPnztW0adPqm17ZNq/i4mJdfvnlrrmXlcuavI0dO9anmZux+9nFnWUIBg8e7J7r6aefbrEeXnzxRe28887q1KmTa0ZmTb2+/fZbn2NOOeUUt2/BggUaPXq0+3/v3r3dRWpFRYXPsZWVlbrlllu05ZZbuufPzs52TdUKCgp8jvM2TXvttde0zTbbuGZs48aNc/usLkaNGuWan9n9zz//fL3zzjs+9W0ZIQtslixZ0uQ1nXbaaeratasLcJuzofekYRkt8zJ8+HBXxs0220z33Xdfk8dr7XtUW1vrmtKNGDHCBbTegObNN9/c4LngPeeeeeYZV/f2HFZ2e1+8Tfca82bC7Lxs/Losm2N1b+Ww88Vue+9jt+2c2HHHHV1TwbayAC8zM9PVmb/ydenSxb3XLbE6GjRokH777Te1F8twWXn8ZbHeffddt8/7nti5e9ZZZ7lz3ntOW7Zx6tSpAf2O2RRbb721rrnmGi1evNh9xhtavny5Cy4POeQQl120czKSMscA2iY0fyoEgI1QU1Ojjz76yP0CbhdugWaBkT2n9cGxplSFhYX6/PPPtXbtWrffgoijjjrK9VOypnrGe1G4fv1692u7/dJ+9dVXu2DDgoDrr79ec+bMcReaDS+eLWNg/Z9sv/W5ycnJafEi1JrKWRBmawvQ7rjjDu2+++766quvtNVWW/lkjOyC1DIyFmBMnz7dBVpWZnsuYxeLhx12mHt+y9BYnyu7MLfmetaf0YKFhpm6mTNn6scff3TPbcGRBRUrVqxwr9f+b83zrPwvvPCCC14bOvvss/XPf/5TDz/8sAsqvayPkzUFtOMtKGvre+JlQZ8FaxZEWX0+99xzuvjii11dWVC3se+RBczPPvusq0cLgOxC3+rBG4S1dC54XXXVVS4wt2A+Nja2xfe4pSDHHscCAXsuC7CPOOIIt+0///lPfYDyj3/8wwWECxcurH/vrKz2flkTypaCA3svv//+ex177LEuYG8rO/fsPLLgqj2DIAtwn3zyySZNJu01WZ3aOWJOPPFE9x7Z+TZw4EB3jtjtVatWhcx3jD/2ebXPoX1WTzrpJJ/XZ+WzH0O8WdInnnjCnQv+AnEAUcIDAGFu5cqVHvs6O+6441p9n5EjR7rF6+OPP3aPYeuGFi5c6LY/+eST7nZhYaG7PX78+BYff8iQIT6P73Xbbbd5YmNjPV9//bXP9ldeecU97pQpU+q32e309HTP6tWrN/h6Fi9e7ImPj/dceOGFPttLSko8PXr08BxzzDH1204++WT32C+99JLPsaNHj/YMGjSo/vYLL7zgjnv11Vd9jrOy2/YJEybUb+vbt68nLi7OM2/ePJ9j//73v3tiYmI8c+fO9dm+//77N6lvK1dOTo6noqKiftu///1vV1/2PjSnte+JldHKMmvWLJ/t++23nyctLc2zbt26jXqPpk+f7m5fc801bToXvOfcHnvs0WTfDTfc4PY1ZuehbW9YH/a6UlJSPEuXLq3fZq/RjuvZs2f96zKvv/662/7mm2/Wb1u0aJF770477bQWX8eXX37p7nvllVd6WsvKZudVVVWVW6zc3vPPzo2Gxx100EF+H8N7vnk/g82577773HENz0H77CQlJXkuu+yy+m2dO3f2jB071hPo75iNfU3e97ygoMDvfcrKytz+Aw88sH5bbW2tZ8CAAZ68vDxPdXW1z+P85z//2YhXCCDSRHVzTfs1zJo3WPMe+7XLfjHfWHYdZh2j7ddA+3XWfuGzX0wBRKasrCxtvvnmLkNmfYqsKaRlvFrLmtANHTrUNe+rrq6uX6xZpb/monvvvbdrHrch77//vnsc+4W/4eNa9suyUo0f157Lvv8asoxVwyZ0VlZrWmfHNXxMK7tlwRo/pt3fvgsbsmaK9nobZhHNX//61yavwTJq+fn5rt+XsXq17N9BBx3k09x1U94TG7zCsj4N/e1vf3PNMy2bszHvkTUDNNb8dFMceeSR2lRWVmvu6WXNM41lXBtm3bzbG77Plvmx1xeoATumTJmihIQEt1jG8KWXXtKFF17ok7FtD8cff7z7O9wwG2lZY2uCbNltL2uyasfY83/55ZcusxgO6n73UZPPlzXvtSxsXFyc22av1c5Ty+YBiF5RHeRZ3wr7Y3///fe3+THsouSxxx5zgd5PP/2kt956y/0BAdBxbJACu5C1JmiBZhdP1vzNLvitaaANnW7Nzi666KIWhzj3shHwbAh270Wvd7G+THYRZ80MG+rZs2eryuUdWW+HHXZo8tjWh6fx41p9NW7+aBfIDfu92WNaUzZrgtj4MVeuXNmqsloTOOsf15i/bdbczpqWPvDAA/XBljUlbNy0c1PeEwtOG/Nu8zbXa+17ZH277MLa32NujNa+xy2xQLche89a2t5S/8bmWDNYs7Gfs912282NRmnNe3/44Qd3TlmfPm9ZjPXHtCaH/lgAauw9aIm9VmvSaNMzeB/Lgjn7m2zBvZd9Hiwosr/d1kzW7mc/jtg53Z7fMe3xmhryBub2w7SXNzA//PDDXb3aYs11rc5fffXVJs2VAUSPqO6TZyOB2dIc66NhfUusz4Z9Udovu//+97/r59ayvif2K7P1UbBO5ACCwy60bZANy6y0dVQ5b8DTeOCRxoGMN/Phvbiyke4sM2F9vOw7w/pVtcQuFq0vVHO/sjccVc+0tk+N9342R5aVrz3YY9qAJ82NGGhBz4bKavf3BqANNXdBbYGZDUdvWTX7Ac4ygzZy4Ya09j3x97zebVbWjXmPLJC0i3i7/6YEav7qreH52LAPn7/zsaPYaxw2bJgbvdH6Lba2X54FHdtvv32Lx1jQb1My+OPd7u+HgcYsi2WZYJvWwYJSCy4bT9Vg79/48ePdYgOZ2IAsV155pcsiN3eut+U7pr1ek5d34BjvNUhRUZEL5Lw/7vjz/PPP67zzzmv1cwCIHFGdyWvNH4v//ve/rtO//aprFx4HHHCA5s+f7/Zb1s5GZrNfm60JijUnOuOMM3wmQwXQMWyACcuynHnmme7CvjFrkmWf2eZ4mwPaZ93fhVVzLAixH4Ps4tfb3M/YhblNIdCYDXrxyy+/uIDCLnwbLy01S2yJZbEsc2CP7e9xN3SR7Y+V1bJbFsj4e7zW/LhlTUXthzDL4DRk36v+WEbCLs5tMBgb4MQuUDd28Ijm3hNjA6g0nmvMLoQtYLUM4Ma8R94fCTc031tz50JLmjsfWzqHO8J1112nNWvWuGDcX/PB0tJSFwRuLBswxN95Yixg906iviE2sqc1W7UBWGyxYNlf02AvO9csU2w/JDQ+Vzb1O6a9XpOxc9a6gth5YZPQe89bO69swKSPP/64yWLBLE02gegV1Zm8ltgfeGvLb7/YeZtG2Mhr9iuf/eGwL1sbstuaT9ivht7mIZdccokbSc1G4QLQcazZlV1sW1BgI+DZfFbWRMsuvKyPlg17btn4xv3QvKzJnV2U3Xbbba4PnGWGrAmgTQnQkF1020Wh/eizxRZbuCZn9nm37ZYN8LIAwwIZaxpmPwbZxaZts5Ed7dd3m0zavi+sH5v1H7OMgl0cW3DT2gu/huziz0Z3tBH17LvJfpCy12FZNBtZ00a39E5p0Fo2JYG1ZLBRCa1pujV7s+Zl9r1oF5E28qYFZS2x12sXmhYQWfksc2EXp9a83dhoko0zJtbHzUaBtDLb6JUb0tr3xNj3uTXpsyyfZaZsZEzL+lgrDW9mqrXvkTUttZEarW+X1bMFhxbQ2flmj2X9zlo6F1pidW7NCL2jdloAb00P/U0xsaLine truncated
                        "text/plain": [
                            "<Figure size 900x600 with 2 Axes>"
                        ]