Refactor ClusterFinderCUDA

Rework the multi-stream pipeline to eliminate per-frame sync barriers and fix the D2H staging architecture. Sync reduction: - Replace one cudaStreamSynchronize per frame with one per stream per batch, cutting synchronisation calls from O(n_frames x n_streams) to O(n_streams) - Introduce a unified per-frame D2H output layout [uint32_t count | clusters[max]] stored in a single class-level lazy-allocated pinned pool (h_output_pinned), replacing the per-stream separate cluster/count device buffers - Move CUDA event pool from per-stream fixed-size to per-frame-slot lazy-allocated, enabling correct kernel timing across any batch size Pinned H2D without CPU-side copy: - Add register_input_buffer(ptr, bytes) / unregister_input_buffer() wrapping cudaHostRegister so callers can pin their existing batch buffer once; all find_clusters_batched() slices then transfer at DMA speed (~22 GB/s) instead of ~15 GB/s for pageable, with no extra memcpy or WC-memory penalty Result (RTX 4090, 400x400 uint16, 3x3 clusters, batch=2000, 5 streams): Before: ~34 µs/frame -> After: ~28 µs/frame (−18 %)
2026-11-03 13:42:28 +01:00 · 2026-05-18 16:30:13 +02:00
1 parent 5cfc888567
commit 6a12e3de24
@@ -55,25 +55,49 @@ def _cuda_available():


 def ClusterFinderCUDA(image_size, cluster_size=(3,3), n_sigma=5, dtype=np.int32,
-                      capacity=1024, n_streams=1):
+                      max_clusters_per_frame=2048, n_streams=4):
    """
    Factory function to create a ClusterFinderCUDA object. Provides a cleaner
    syntax for the templated ClusterFinderCUDA in C++. API mirrors
-    ClusterFinder() plus CUDA-specific knobs (n_streams).
+    ClusterFinder() plus CUDA-specific knobs.

+    Parameters
+    ----------
+    image_size : tuple of (int, int)
+        Detector shape as (nrows, ncols).
+    cluster_size : tuple of (int, int), optional
+        Cluster window size; default (3, 3).
+    n_sigma : float, optional
+        Threshold in units of per-pixel pedestal standard deviation.
+    dtype : numpy dtype, optional
+        Cluster value type (np.int32 or np.float32).
+    max_clusters_per_frame : int, optional
+        Tight upper bound on clusters per frame. Determines the fixed-size D2H
+        transfer per frame. Set this high enough to never truncate real frames
+        but as tight as possible to minimize PCIe traffic. Default 2048.
+    n_streams : int, optional
+        Number of CUDA streams for H2D/kernel/D2H pipelining. Default 4.
+
+    Example
+    -------
    .. code-block:: python

        from aare import ClusterFinderCUDA

-        cf = ClusterFinderCUDA(image_size=(512, 1024),
+        cf = ClusterFinderCUDA(image_size=(400, 400),
                               cluster_size=(3, 3),
                               n_sigma=5,
-                               n_streams=4)
+                               n_streams=5)
        for frame in pedestal_frames:
            cf.push_pedestal_frame(frame)
+
+        # Batched (recommended for throughput)
+        results = cf.find_clusters_batched(frames_3d, first_frame=0)
+
+        # Or single-frame (one launch per frame)
        for i, frame in enumerate(data_frames):
            cf.find_clusters(frame, frame_number=i)
-        clusters = cf.steal_clusters()
+            clusters = cf.steal_clusters()
    """
    if not _cuda_available():
        raise RuntimeError(
@@ -84,7 +108,7 @@ def ClusterFinderCUDA(image_size, cluster_size=(3,3), n_sigma=5, dtype=np.int32,
    cls = _get_class("ClusterFinderCUDA", cluster_size, dtype)
    return cls(image_size,
               n_sigma=n_sigma,
-               capacity=capacity,
+               max_clusters_per_frame=max_clusters_per_frame,
               n_streams=n_streams)

 def ClusterCollector(clusterfindermt, dtype=np.int32): 
@@ -8,8 +8,8 @@

 #include <cstdint>
 #include <pybind11/pybind11.h>
-// #include <pybind11/stl.h>
-#include <pybind11/stl_bind.h>
+#include <pybind11/stl.h>
+// #include <pybind11/stl_bind.h>

 namespace py = pybind11;
 using pd_type = double;
@@ -28,18 +28,20 @@ void define_ClusterFinderCUDA(py::module &m, const std::string &typestr) {

    using ClusterType = Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>;
    using CF = ClusterFinderCUDA<ClusterType, uint16_t, pd_type>;
+    using ContigArr =
+        py::array_t<uint16_t, py::array::c_style | py::array::forcecast>;

    py::class_<CF>(m, class_name.c_str())
-        .def(py::init<Shape<2>, pd_type, size_t, int>(), py::arg("image_size"),
-             py::arg("n_sigma") = 5.0, py::arg("capacity") = 1'000'000,
-             py::arg("n_streams") = 1)
+        .def(py::init<Shape<2>, float, size_t, int>(), py::arg("image_size"),
+             py::arg("n_sigma") = 5.0f,
+             py::arg("max_clusters_per_frame") = 2048, py::arg("n_streams") = 4)

        .def_property(
            "nSigma", &CF::get_nSigma, &CF::set_nSigma,
            R"(Number of sigma above the pedestal to consider a photon during cluster finding.)")

        .def("push_pedestal_frame",
-             [](CF &self, py::array_t<uint16_t> frame) {
+             [](CF &self, ContigArr frame) {
                 auto view = make_view_2d(frame);
                 self.push_pedestal_frame(view);
             })
@@ -63,39 +65,59 @@ void define_ClusterFinderCUDA(py::module &m, const std::string &typestr) {
        .def(
            "steal_clusters",
            [](CF &self, bool realloc_same_capacity) {
-                ClusterVector<ClusterType> clusters =
-                    self.steal_clusters(realloc_same_capacity);
-                return clusters;
+                return std::move(self.steal_clusters(realloc_same_capacity));
            },
-            py::arg("realloc_same_capacity") = false)
+            py::arg("realloc_same_capacity") = true)

        .def(
            "find_clusters",
-            [](CF &self, py::array_t<uint16_t> frame, uint64_t frame_number) {
+            [](CF &self, ContigArr frame, uint64_t frame_number) {
                auto view = make_view_2d(frame);
                self.find_clusters(view, frame_number);
            },
-            py::arg("frame"), py::arg("frame_number") = 0)
+            py::arg("frame"), py::arg("frame_number") = 0,
+            py::call_guard<py::gil_scoped_release>())

        .def(
            "find_clusters_batched",
-            [](CF &self, py::array_t<uint16_t> frames, uint64_t first_frame) {
-                // frames is expected as a 3D numpy array (n_frames, nrows,
-                // ncols)
+            [](CF &self, ContigArr frames, uint64_t first_frame) {
                auto view = make_view_3d(frames);
                return self.find_clusters_batched(view, first_frame);
            },
            py::arg("frames"), py::arg("first_frame") = 0,
-            R"(Process a 3D array of frames (n_frames, nrows, ncols) in parallel
-            across the configured CUDA streams. Returns a list of ClusterVector, one per
-            input frame.)")
+            py::call_guard<py::gil_scoped_release>(),
+            R"(Process a 3D array of frames (n_frames, nrows, ncols) using
+            n_streams CUDA streams for H2D/kernel/D2H pipelining. Returns a
+            list of ClusterVector, one per input frame. The input array is
+            converted to C-contiguous uint16 if needed.)")

        .def("avg_kernel_time_ms", &CF::avg_kernel_time_ms,
             R"(Average kernel execution time per frame in milliseconds,
-            excluding PCIe transfers. Use wall_time - avg_kernel_time to estimate transfer overhead.)")
+            excluding PCIe transfers.)")

        .def("reset_timers", &CF::reset_timers,
-             R"(Reset the internal kernel timing counters.)");
+             R"(Reset the internal kernel timing counters.)")
+
+        .def(
+            "register_input_buffer",
+            [](CF &self, py::array arr) {
+                auto info = arr.request();
+                self.register_input_buffer(
+                    info.ptr, static_cast<size_t>(info.size) *
+                                  static_cast<size_t>(info.itemsize));
+            },
+            R"(Pin a numpy array as a locked host buffer so that
+            find_clusters_batched transfers it at full DMA bandwidth
+            (~22 GB/s) instead of going through the CUDA driver's
+            internal staging (~15 GB/s for pageable memory).
+
+            Call once before the processing loop with the full data
+            array.  Slices of that array passed to find_clusters_batched
+            lie within the registered region and benefit automatically.
+            Call unregister_input_buffer() when done.)")
+
+        .def("unregister_input_buffer", &CF::unregister_input_buffer,
+             "Release the previously pinned input buffer.");
 }

 } // namespace aare
@@ -53,7 +53,7 @@
                    "text": [
                        "Image size:       (400, 400)\n",
                        "Pedestal frames:  1000\n",
-                        "Data frames:      88999\n"
+                        "Data frames:      20000\n"
                    ]
                }
            ],
@@ -62,10 +62,10 @@
                "f = File(base / 'Moench03new/cu_half_speed_master_4.json')\n",
                "\n",
                "n_frames_pd = 1000\n",
-                "N           = 88999\n",
+                "N           = 20000 #88999\n",
                "cluster_size = (9, 9)\n",
                "image_size   = (f.rows, f.cols)\n",
-                "capacity     = 1_000 #3_000_000\n",
+                "capacity     = 50_000 #3_000_000\n",
                "\n",
                "print(f'Image size:       {image_size}')\n",
                "print(f'Pedestal frames:  {n_frames_pd}')\n",
@@ -136,17 +136,33 @@
        {
            "cell_type": "code",
            "execution_count": 7,
-            "id": "cbdcb805-708b-4205-bda9-2aa163d0e81f",
+            "id": "00e2f63a-5e64-4d3e-becf-5e2835c6d712",
            "metadata": {},
            "outputs": [],
            "source": [
-                "N_STREAMS  = 5\n",
-                "cf_cuda = ClusterFinderCUDA(image_size, cluster_size, n_sigma=7, capacity=capacity, n_streams=N_STREAMS)"
+                "# Runs the destructor under the hood in case cf_cuda has already been constructed\n",
+                "# del cf_cuda\n",
+                "cf_cuda = None"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 8,
+            "id": "cbdcb805-708b-4205-bda9-2aa163d0e81f",
+            "metadata": {},
+            "outputs": [],
+            "source": [
+                "N_STREAMS  = 5\n",
+                "cf_cuda = ClusterFinderCUDA(image_size, \n",
+                "                            cluster_size, \n",
+                "                            n_sigma=7, \n",
+                "                            max_clusters_per_frame=2048,\n",
+                "                            n_streams=N_STREAMS) "
+            ]
+        },
+        {
+            "cell_type": "code",
+            "execution_count": 9,
            "id": "1546f405-1bf6-4073-ab35-8134be695a6c",
            "metadata": {},
            "outputs": [
@@ -154,7 +170,7 @@
                    "name": "stdout",
                    "output_type": "stream",
                    "text": [
-                        "Pedestal (1000 frames): 0.515s\n"
+                        "Pedestal (1000 frames): 0.496s\n"
                    ]
                }
            ],
@@ -177,7 +193,7 @@
        },
        {
            "cell_type": "code",
-            "execution_count": 9,
+            "execution_count": 10,
            "id": "4573be3d-5ba8-4e18-bab0-c874f2b7dcb2",
            "metadata": {},
            "outputs": [
@@ -185,7 +201,7 @@
                    "name": "stdout",
                    "output_type": "stream",
                    "text": [
-                        "Reading 88999 frames:        4.574s  (19458 FPS, 5938.097 GB/s)\n"
+                        "Reading 20000 frames:        1.195s  (16742 FPS, 5109.275 GB/s)\n"
                    ]
                }
            ],
@@ -208,7 +224,7 @@
        },
        {
            "cell_type": "code",
-            "execution_count": 10,
+            "execution_count": 11,
            "id": "fbb14fda-2852-4b73-ba2c-9952abac1d99",
            "metadata": {},
            "outputs": [
@@ -216,7 +232,7 @@
                    "name": "stdout",
                    "output_type": "stream",
                    "text": [
-                        "CPU clustering:          837.700s (106 FPS, 90799856 clusters, 1020.23/frame)\n"
+                        "CPU clustering:          188.615s (106 FPS, 20325971 clusters, 1016.30/frame)\n"
                    ]
                }
            ],
@@ -257,7 +273,7 @@
        },
        {
            "cell_type": "code",
-            "execution_count": 19,
+            "execution_count": 12,
            "id": "4f94cf35-5796-463d-bed8-f44a02d91fc6",
            "metadata": {},
            "outputs": [],
@@ -267,7 +283,7 @@
        },
        {
            "cell_type": "code",
-            "execution_count": 20,
+            "execution_count": 45,
            "id": "50848f8c-2e66-45b8-b6d5-ca80f368e6ba",
            "metadata": {},
            "outputs": [],
@@ -275,6 +291,9 @@
                "if(BATCHED):\n",
                "    BATCH_SIZE = 5000\n",
                "\n",
+                "    # Before warmup, pin the entire data array once\n",
+                "    cf_cuda.register_input_buffer(data)\n",
+                "\n",
                "    # Warmup: first kernel launch pays CUDA context + pedestal H2D upload cost\n",
                "    _ = cf_cuda.find_clusters_batched(data[0:BATCH_SIZE], first_frame=0)\n",
                "    \n",
@@ -288,14 +307,17 @@
                "            cf_cuda.find_clusters_batched(data[start:stop], first_frame=start)\n",
                "        )\n",
                "    t_cuda = time.perf_counter() - t0\n",
+                "\n",
+                "    cf_cuda.unregister_input_buffer()  # release when done with this dataset\n",
+                "\n",
                "    kernel_ms = cf_cuda.avg_kernel_time_ms()\n",
                "    \n",
                "    n_clusters_cuda = sum(cv.size for cv in clusters_cuda_per_frame)\n",
                "\n",
                "    hist_cuda = make_hist_from_batch(clusters_cuda_per_frame)\n",
                "    \n",
-                "else:\n",
-                "    # Warmup\n",
+                "else:   \n",
+                "    # Simpler: (non-batched) per-frame run on non-pinned data\n",
                "    cf_cuda.find_clusters(data[0])\n",
                "    _ = cf_cuda.steal_clusters(realloc_same_capacity=False)\n",
                "\n",
@@ -322,7 +344,7 @@
        },
        {
            "cell_type": "code",
-            "execution_count": 21,
+            "execution_count": 46,
            "id": "502d0d3b-6b1e-4cc9-91df-9d998bd849b5",
            "metadata": {},
            "outputs": [
@@ -332,7 +354,7 @@
                            "(9, 9)"
                        ]
                    },
-                    "execution_count": 21,
+                    "execution_count": 46,
                    "metadata": {},
                    "output_type": "execute_result"
                }
@@ -343,7 +365,7 @@
        },
        {
            "cell_type": "code",
-            "execution_count": 22,
+            "execution_count": 47,
            "id": "2e3e4b9c-7f23-4fb7-bc3e-fa3d766d51fc",
            "metadata": {},
            "outputs": [
@@ -351,7 +373,7 @@
                    "name": "stdout",
                    "output_type": "stream",
                    "text": [
-                        "CPU clustering:          837.700s (106 FPS, 90799856 clusters, 1020.23/frame)\n"
+                        "CPU clustering:          188.615s (106 FPS, 20325971 clusters, 1016.30/frame)\n"
                    ]
                }
            ],
@@ -362,7 +384,7 @@
        },
        {
            "cell_type": "code",
-            "execution_count": 23,
+            "execution_count": 48,
            "id": "4b8df93b-9a1b-41a5-9fed-9cda295fb523",
            "metadata": {},
            "outputs": [
@@ -370,10 +392,10 @@
                    "name": "stdout",
                    "output_type": "stream",
                    "text": [
-                        "CUDA clustering:          6.261s  (14214 FPS, 88537254 clusters, 994.81/frame)\n",
-                        "  Kernel only:            0.040 ms/frame\n",
-                        "  PCIe + overhead:        0.030 ms/frame\n",
-                        "Speedup (CPU / CUDA):     133.79×\n"
+                        "CUDA clustering:          1.304s  (15332 FPS, 20238568 clusters, 1011.93/frame)\n",
+                        "  Kernel only:            0.038 ms/frame\n",
+                        "  PCIe + overhead:        0.027 ms/frame\n",
+                        "Speedup (CPU / CUDA):     144.60×\n"
                    ]
                }
            ],
@@ -397,7 +419,7 @@
        },
        {
            "cell_type": "code",
-            "execution_count": 24,
+            "execution_count": 49,
            "id": "d3a850df-7df0-485b-971d-381cfdc6be81",
            "metadata": {},
            "outputs": [
@@ -405,7 +427,7 @@
                    "name": "stdout",
                    "output_type": "stream",
                    "text": [
-                        "Cluster count diff:       2262602 (2.49%)\n"
+                        "Cluster count diff:       87403 (0.43%)\n"
                    ]
                }
            ],
@@ -425,7 +447,7 @@
        },
        {
            "cell_type": "code",
-            "execution_count": 25,
+            "execution_count": 50,
            "id": "aff8db7e-8332-4228-857f-a9875cf940d3",
            "metadata": {},
            "outputs": [
@@ -445,14 +467,14 @@
        },
        {
            "cell_type": "code",
-            "execution_count": 26,
+            "execution_count": 51,
            "id": "9adeea2f-9309-4c0a-905b-7d1203ba1f4e",
            "metadata": {},
            "outputs": [
                {
                    "data": {
                        "image/png":
-                            "iVBORw0KGgoAAAANSUhEUgAAAxYAAAJOCAYAAAAqFJGJAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAr9JJREFUeJzs3Qd4VMXaB/B/ek9IIaTQewcFRbAAgiAKiA0rgiKg9AtYEJVyVfwsgA2wgmIBvZSrgAgIglxAkCJFqtJbEtJ7O9/zznqW3c2mhy3Z/+95ls3uzp5zdnI2zHvmnRk3TdM0EBERERERVYJ7Zd5MRERERETEwIKIiIiIiKoEeyyIiIiIiKjSGFgQEREREVGlMbAgIiIiIqJKY2BBRERERESVxsCCiIiIiIgqjYEFERERERFVGgMLIiIiIiKqNAYWRFQl9u3bh8cffxwNGjSAr68vAgMDce211+KNN95AYmKisVy3bt3U7WqZO3cuFi5ceNW2T86jOp4Lv/76KwYOHIjY2Fh4e3sjJCQEXbp0wbx585CRkWEsV79+fbi5uRlv8n3s1KkTvvjiC7PtSbm+ffta3dfvv/+u3usodVjWvzHl/UzTpk0zqyt/f3/Url0bvXv3xnvvvYe0tLRijykvLw9RUVHqff/5z3+q+BMTOR8GFkRUaR9//DE6dOiAnTt34plnnsGaNWuwfPly3H///Zg/fz6GDh1qs1qujo1Jqpjqdi5MnToVt9xyC86dO4d///vfWLduHRYvXowePXqoxvGLL75oVv7GG2/Etm3b1E3qQRq/gwcPVkGIs7HF3xjZptSV3L/11luoW7cunn32WbRq1Qp//PGH1fesXLkSly5dUj9/+umnlT4GImfnae8DICLnJv8RP/3007jtttuwYsUK+Pj4GF+T5yZOnKj+o3ZmmqYhOzsbfn5+cAZZWVlOc6yOQq48S8Pb09Mx/1v87rvvMGPGDNWAlka2HKuuT58+qgEs30VTNWrUwA033GB83LNnT9SrVw+zZs1S31lnYau/MRK4REREGB8/+OCDGD16NLp27Yr+/fvj6NGjZvvWgwnpOZIya9euxdmzZ1VvB5GrYo8FEVXKa6+9pho5H330UZH/dIX8pyv/KRfnl19+Ue+Xe1MnT54skrLw999/q//sY2Ji1L5q1aqlrtbu3bvXmAJx8OBBbNq0yZjWIM/pUlNTMWnSJJVKIccl6STjx483SyER8j5pUMiV0BYtWqh9ff755yXWw5IlS9C5c2cEBASoFA1Jo9izZ49ZmSFDhqjXjh8/jjvuuEP9XKdOHdUwysnJMSubm5uLV155Bc2bN1f7r1mzpkoDiY+PNyunp30sW7YM11xzjUoRmT59unpN6qJXr14qtUPeP2rUKKxatcqsvuXKtzSmz5w5U+QzPfHEEwgPD1dBVXFK+52YHqNcYW7btq06xoYNG+Ldd98tsr2y/o4KCwtVmkr79u1VEKU3or///vtSzwX9nFu0aJGqe9mHHLv8XvS0GEv6FX85Ly0/l1y1lrqX45DzRR7r75HHck5cf/31Kg2noiSoCA0NVXVm7fiCgoLU77okUkfNmjXDqVOnUFXkSr4cj7Wr9T/++KN6Tf+dyLk7fPhwdc7r57T0qqxfv/6q/o2pjHbt2mHKlCk4ffq0+o6bOn/+vApo+vXrp3pR5JysTj1kRBXhmJdmiMgpFBQUYMOGDepKnzQWrjZpjMs+Jada0hQSEhKwdetWJCcnq9el4XrfffepvHNJgxF6QyQzM1NdVZQrii+88IJq4ErD8+WXX8b+/ftV48a0wSZXRiWfXV6XHOrIyMgSGz6ShiINf7mXoODNN9/EzTffjB07dqBly5ZmV8alESRXnqVRu3nzZtW4l2OWfQlpoNx1111q/3IlWnLopTEoqTAyPkUaqKY9Ert378ahQ4fUvqVBLg3ZCxcuqM8rP0vqixz/N998owImUyNGjMCrr76KDz/8UAUyOslZlzQbKS+BQEV/JzoJNCRAkIa71OdXX32FcePGqbqSQKK8vyMJ0r788ktVj9Lolsal1IPe8C/pXNBNnjxZBYMSQLq7u5f4Oy6pYS3bkcan7EuCunvuuUc99/PPPxsbxc8995wKQk6cOGH83cmxyu9L0pNKapDK7/LAgQN44IEHVJBYUXLuyXkkDfqqbHhLULVgwYIi6UjymaRO5RwRgwYNUr8jOd+aNm2qzhF5fPnyZYf5G2ONfF/leyjf1ccee8zs88nxSQCu9wZ99tln6lywFvwRuQSNiKiCLl68qMmfkQcffLDM7+natau66TZu3Ki2IfemTpw4oZ5fsGCBepyQkKAez5kzp8Ttt2rVymz7upkzZ2ru7u7azp07zZ7/z3/+o7a7evVq43PyOCQkREtMTCz185w+fVrz9PTUxowZY/Z8WlqaFhUVpQ0cOND43ODBg9W2v/32W7Oyd9xxh9asWTPj42+++UaVW7p0qVk5OXZ5fu7cucbn6tWrp3l4eGhHjhwxK/vMM89obm5u2sGDB82e7927d5H6luOKjIzUcnJyjM/93//9n6ov+T0Up6y/EzlGOZa9e/eaPX/bbbdpwcHBWkZGRrl+R5s3b1aPp0yZUqFzQT/nbrnlliKvTZ06Vb1mSc5Ded60PuRz+fn5aWfPnjU+J59RykVHRxs/l1ixYoV6/vvvvzc+d/LkSfW7e+KJJ0r8HNu3b1fvff7557WykmOT8yovL0/d5Lj180/ODdNyd955p9Vt6Oeb/h0szrvvvqvKmZ6D8t3x8fHRJk6caHwuMDBQGz9+vHa1/8aU9zPpv/P4+Hir78nKylKv9+nTx/hcYWGh1rhxYy02NlbLz883287PP/9cjk9IVL0wFYqInEJYWBgaNWqkegIkR1zSjOTKfllJekrr1q1V6kx+fr7xJilL1lKxbr31VpV6UpqffvpJbUeuZJpuV67yy9V3y+3KviR1wpRcmTdNT5FjlbQVKWe6TTl2udpvuU15v1wBNiUpQPJ5TXtLxEMPPVTkM0jPQVxcnMrjF1Kv0stx5513mqWSVeZ3IgNg5eq2qYcfflilPslV6/L8jiTFRkhqV2Xce++9qCw5Vkml0knqk5CeJdPeBf1509+zXOGWz3e1Bv2uXr0aXl5e6iY9I99++y3GjBlj1jNVFR555BHVG2Ta6yK9Y5LeJ714OkkHkzKy/+3bt6seFGdguNaAIt8vSZ2T3iYPDw/1nHxWOU+l14LIVTGwIKIKk4GO0niS9I6rTf7DltQSaWRK2o1MMykpHWPHji1xOkidzNwi01XqDS39Jrnp0nCQFB5T0dHRZToufUaY6667rsi2JSfbcrtSX5apRdIoMx3HINuUNBFJ77Hc5sWLF8t0rJJeIuMdLFl7TlJZJG3rgw8+MDbwJU3HMm2qMr8TCYgs6c/pqTBl/R1Jrr405qxtszzK+jsuiQRXpuR3VtLzJY1XKY6kmInyfs9uuukmNYuSpM79+eef6pySMRr6sQgZXyPpPNZI0CPkd1AS+aySLiRT2erbkgBCAgkJKHXyfZCG+CeffKJS0OR9EpDLOV2Vf2Oq4jOZ0oNBGUek04PBu+++W9Wr3CQVTup86dKlRVIBiVwFx1gQUYVJ404G6soV5IrOhqI3si0HL1s2nvUrvPp/6DJDi1yBlZx9ydOXPPmSSANFctuLu5poOhuMKGuOtP4+mcNejq8qyDZl0HRxM91IQ7u0Y5X360GPqeIacRIMyNSd0nvw/vvvqx4QmXGnNGX9nVjbr/6cHGt5fkcSvEjDUd5fmeDAWr2Zno+mYzKsnY+2Ip+xTZs2atYhGYdS1nEW0tDt2LFjiWUk0JTpa63Rn7cWjFqSq/XS4yVT4EogJAGN5bS28vubM2eOuslgaBnU/fzzz6vesuLO9Yr8jamqz6TTB5/r6++kpKSo4EG/oGDN119/jZEjR5Z5H0TVBXssiKhSZJCqXE0eNmyYakxaknSHH374odj366k2cqXa2n/mxZGGrwxWlgaXnkojpDEo061akoGzf/31l2rESmPL8lZSyk9J5Gq9XCGVbVvbbmkNO2vkWOUqvjSerW1PZvYpjaRhyYBfuVJtSgZkWyNXXqVBKAPKZZC0NIrKOwC1uN+JkEHYlmsBSONLgiTp6SjP70imVxWlrcdQ3LlQkuLOx5LOYVt46aWXkJSUpAJAa6k56enpKvAoLxl0bO08ERIk6gvrlUZmpJKUMBnELTcJ0Kyl3enkXJMeMQleLc+Vyv6NqarPJOSclQH4cl7IwoT6eSvnlUy6sHHjxiI3CaCYDkWuij0WRFQpktIgDTxpiMrMLTLfvKQ/yH/2knMvU0RK3rzluAKdpLNIQ2DmzJlqTINcAZf0Gpk+1ZQ09KQhIlfVmzRpotI5ZLYYeV6ueuqkUSuNZ0m7kClNpYEjz8mMRHKVURYY+9e//qXGJch4ALlyKg0yaVCXtbFhShocMiuRzAQjU6/efvvt6nNIb4HMCCWzMunTv5aVTN8qsybJbDoy/kFSSiR1Q67YSsNFZoySQKAk8nmlcSONcDk+uUIrDaLDhw+r12UWJMsrwzJmQWYvkmOWWZdKU9bfiZ5GIuky0pshV+BlRie5uv1///d/xivwZf0dSdqWzDAkufpSzxKQSBAh55tsS8YRlHQulETqXFJ09NmmJGiUtB5r0/FWlqTYyBgVSQ8qbZyF1LEEF9KYld+hHJ+8V3owfvvtNzWrl8waVdqUs5bk/JIUJrkaLzNxSf1IACN1Jr1wMnbGsofMGjl/JK1JygcHB6uZsaTHRCdX+bt3767G1cgUyrJN6dWQngopW5V/Yyr6mXbt2qWOWbYrU8nK3yGZklhmtpLARU8hk9+VfMdlNjNrM6bp9SBBieW4IqJqz96jx4moepDZcGTWmbp162re3t5aQECAds0112gvv/yyFhcXV+ysUOLChQvafffdp4WFhanZmB599FHt999/N5u95dKlS9Line truncated
+                            "iVBORw0KGgoAAAANSUhEUgAAAxYAAAJOCAYAAAAqFJGJAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAtShJREFUeJzs3Qd4FNX6BvA3IY1ACITeRXpTlCIgCgqCCja8YgNREVQQRcSC2PCqWBCwoIgNL6ioF/GiIIIN5Q8IUkSkikgvAUJoIQlk/s971lk2m03Z7CbZ8v6eZ9lkdnZ29uxsmG/Od74TYVmWBRERERERER9E+vJkERERERERBRYiIiIiIuIX6rEQERERERGfKbAQERERERGfKbAQERERERGfKbAQERERERGfKbAQERERERGfKbAQERERERGfKbAQERERERGfKbAQEb9YvXo1brvtNtSrVw9xcXEoW7Yszj33XLz44os4ePCgc70uXbqYW1F54403MGXKlCLbvgSPUDwWfv75Z/Tp0wc1a9ZETEwMEhMT0bFjR7z55ps4duyYc70zzjgDERERzhu/j+eddx7+85//ZNse1+vVq5fH1/r111/NcwOlDQv6N8bb9/TUU09la6v4+HjUqlULPXr0wGuvvYYjR47kuk+ZmZmoVq2aed5///tfP79jkeCjwEJEfPb222+jdevWWLZsGR588EHMnTsXM2fOxHXXXYdJkyZhwIABxdbKoXgyKYUTasfCk08+iQsvvBA7d+7Ev//9b8yfPx/Tp09H165dzcnxY489lm39888/H4sXLzY3tgNPfvv372+CkGBTHH9juE22Fe/Hjh2LOnXq4KGHHkLz5s3x22+/eXzOV199hb1795qf3333XZ/3QSTYRZX0DohIcON/xHfffTcuueQSfPHFF4iNjXU+xmUPPPCA+Y86mFmWhRMnTqB06dIIBmlpaUGzr4GCV5554h0VFZj/LX722Wd4+umnzQk0T7K5r7bLLrvMnADzu+iqfPnyaN++vfP3bt26oW7duhg3bpz5zgaL4vobw8ClUqVKzt9vuOEG3HPPPejcuTOuvPJKbNy4Mdtr28EEe464zrx587Bjxw7T2yESrtRjISI+ee6558xJzuTJk3P8p0v8T5f/Kefmxx9/NM/nvau///47R8rCX3/9Zf6zr1GjhnmtqlWrmqu1q1atcqZA/PHHH1iwYIEzrYHLbIcPH8aIESNMKgX3i+kkw4YNy5ZCQnweTyh4JbRp06bmtT744IM82+GTTz5Bhw4dUKZMGZOiwTSKlStXZlvn1ltvNY/9+eefuPzyy83PtWvXNidG6enp2dbNyMjAM888gyZNmpjXr1y5skkDSU5Ozraenfbx+eef45xzzjEpIqNHjzaPsS26d+9uUjv4/CFDhmD27NnZ2ptXvnkyvX379hzv6fbbb0fFihVNUJWb/D4T133kFeazzjrL7OOZZ56JV199Ncf2CvoZZWVlmTSVVq1amSDKPomeNWtWvseCfcxNnTrVtD1fg/vOz8VOi3FnX/Hncen+vnjVmm3P/eDxwt/t5/B3HhPt2rUzaTiFxaCiQoUKps087V9CQoL5rPPCNmrcuDG2bt0Kf+GVfO6Pp6v1X3/9tXnM/kx47A4aNMgc8/YxzV6Vb7/9tkj/xvji7LPPxqhRo7Bt2zbzHXe1a9cuE9BcccUVpheFx2Qo9ZCJFEZgXpoRkaBw6tQpfP/99+ZKH08WihpPxvmazKlmmsL+/fuxaNEiHDp0yDzOE9d//etfJu+caTBkn4gcP37cXFXkFcVHH33UnODyxPOJJ57A77//bk5uXE/YeGWU+ex8nDnUVapUyfPEh2koPPHnPYOCl156CRdccAGWLl2KZs2aZbsyzpMgXnnmSe1PP/1kTu65z3wt4gnKVVddZV6fV6KZQ8+TQabCcHwKT1BdeyRWrFiBdevWmdfmCTlPZHfv3m3eL39m6gv3/+OPPzYBk6s777wTzz77LN566y0TyNiYs840G67PQKCwn4mNgQYDBJ64sz0//PBD3HfffaatGEh4+xkxSJs2bZppR5508+SS7WCf+Od1LNhGjhxpgkEGkJGRkXl+xnmdWHM7PPnkazGo6927t1n23XffOU+KH374YROEbNmyxfnZcV/5eTE9Ka8TUn6Wa9aswfXXX2+CxMLiscfjiCf0/jzxZlD1/vvv50hH4ntim/IYoX79+pnPiMdbo0aNzDHC3w8cOBAwf2M84feV30N+V2+55ZZs74/7xwDc7g167733zLHgKfgTCQuWiEgh7dmzx+KfkRtuuKHAz+ncubO52X744QezDd672rJli1n+/vvvm9/3799vfp8wYUKe22/evHm27dvGjBljRUZGWsuWLcu2/L///a/Z7pw5c5zL+HtiYqJ18ODBfN/Ptm3brKioKGvo0KHZlh85csSqVq2a1adPH+ey/v37m21/+umn2da9/PLLrcaNGzt///jjj816M2bMyLYe953L33jjDeeyunXrWqVKlbI2bNiQbd0HH3zQioiIsP74449sy3v06JGjvblfVapUsdLT053LXnjhBdNe/BxyU9DPhPvIfVm1alW25ZdccolVrlw569ixY159Rj/99JP5fdSoUYU6Fuxj7sILL8zx2JNPPmkec8fjkMtd24Pvq3Tp0taOHTucy/geuV716tWd74u++OILs3zWrFnOZX///bf57G6//fY838eSJUvMcx955BGroLhvPK4yMzPNjfttH388NlzX69mzp8dt2Meb/R3MzauvvmrWcz0G+d2JjY21HnjgAeeysmXLWsOGDbOK+m+Mt+/J/syTk5M9PictLc08ftlllzmXZWVlWQ0aNLBq1qxpnTx5Mtt2vvvuOy/eoUhoUSqUiASFpKQk1K9f3/QEMEecaUa8sl9QTE9p0aKFSZ05efKk88aUJU+pWBdffLFJPcnPN998Y7bDK5mu2+VVfl59d98uX4upE654Zd41PYX7yrQVrue6Te47r/a7b5PP5xVgV0wB4vt17S2hG2+8Mcd7YM/Bvn37TB4/sV3Zy9GzZ89sqWS+fCYcAMur265uuukmk/rEq9befEZMsSGmdvni2muvha+4r0ylsjH1idiz5Nq7YC93/Zx5hZvvr6gG/c6ZMwfR0dHmxp6RTz/9FEOHDs3WM+UPN998s+kNcu11Ye8Y0/vYi2djOhjX4esvWbLE9KAEA8e1BuT4fjF1jr1NpUqVMsv4XnmcstdCJFwpsBCRQuNAR548Mb2jqPE/bKaW8CSTaTcsM8mUjnvvvTfPcpA2Vm5huUr7RMu+MTedJw5M4XFVvXr1Au2XXRGmbdu2ObbNnGz37bK93FOLeFLmOo6B22SaCNN73Le5Z8+eAu0r00s43sGdp2VMZWHa1sSJE50n+EzTcU+b8uUzYUDkzl5mp8IU9DNirj5P5jxt0xsF/YzzwuDKFT+zvJbnNV4lN0wxI2+/Z506dTJVlJg6t3btWnNMcYyGvS/E8TVM5/GEQQ/xM8gL3yvThVjK1t4WAwgGEgwobfw+8ET8nXfeMSlofB4Dch7T/vwb44/35MoOBjmOyGYHg9dcc41pV96YCsc2nzFjRo5UQJFwoTEWIlJoPLnjQF1eQS5sNRT7JNt98LL7ybN9hdf+D50VWngFljn7zNNnnnxeeILC3Pbcria6VoOhguZI289jDXvunz9wmxw0nVulG55o57evfL4d9LjK7SSOwQBLd7L34PXXXzc9IKy4k5+CfiaeXtdexn315jNi8MITRz7fl+DAU7u5Ho+uYzI8HY/Fhe+xZcuWpuoQx6EUdJwFT3TbtGmT5zoMNFm+1hN7uadg1B2v1rPHiyVwGQgxoHEva8vPb8KECebGwdAc1P3II4+Y3rLcjvXC/I3x13uy2YPP7fl3UlNTTfBgX1Dw5KOPPsLgwYML/BoioUI9FiLiEw5S5dXkgQMHmpNJd0x3+PLLL3N9vp1qwyvVnv4zzw1PfDlYmSdcdioN8WSQ5VbdceDs5s2bzUksT7bcb3ml/OSFV+t5hZTb9rTd/E7sPOG+8io+T549bY+VffLDNCwO+OWValcckO0Jr7zyhJADyjlImidF3g5Aze0zIQ7Cdp8LgCdfDJLY0+HNZ8TyqpTffAy5HQt5ye14zOsYLg6PP/44UlJSTADoKTXn6NGjJvDwFgcdezpOiEGiPbFefliRiilhHMTNGwM0T2l3Nh5r7BFj8Op+rPj6N8Zf74l4zHIAPo8LTkxoH7c8rlh04YcffshxYwCldCgJV+qxEBGfMKWBJ3g8EWXlFtabZ/oD/7Nnzj1LRDJv3n1cgY3pLDwRGDNmjBnTwCvgTK9h+VRXPNHjiQivqjds2NCkc7BaDJfzqqeNJ7U8eWbaBUua8gSHy1iRiFcZOcHY/fffb8YlcDwAr5zyhIwn1AU92XDFEw5WJWIlGJZevfTSS837YG8BK0KxKpNd/rWgWL6VVZNYTYfjH5hSwtQNXrHliQsrRjEQyAvfL09ueBLO/eMVWp4QrV+/3jzOKkjuV4Y5ZoHVi7jPrLqUn4J+JnYaCdNl2JvBK/Cs6MSr2y+88ILzCnxBPyOmbbHCEHP12c4MSBhE8HjjtjiOIK9jIS9sc6bo2NWmGDQyrcdTOV5fMcWGY1SYHpTfOAu2MYMLnszyM+T+8bnswfjll19MVS9Wjcqv5Kw7Hl9MYeLVeFbiYvswgGGbsReOY2fce8g84fHDtCauX65cOVMZiz0mNl7lv+iii8y4GpZQ5jbZq8GeCq7rz78xhX1Py5cvN/vM7bKULP8OsSQxK1sxcLFTyPhZ8TvOamaeKqbZ7cCgxH1ckUjIK+nR4yISGlgNh1Vn6tSpY8XExFhlypSxzjnnHOuJJ56w9u3bl2tVKNq9e7f1r3/9y0pKSjLVmPr27Wv9+uuv2aq37N2717r11lutJk2amG2zwsxZZ51ljR8/3lmVxa600717dyLine truncated
                        "text/plain": [
                            "<Figure size 800x600 with 2 Axes>"
                        ]
@@ -495,7 +517,7 @@
        {
            "cell_type": "code",
            "execution_count": null,
-            "id": "e33bdd78-8f0c-4690-9153-3f9eaffd4afa",
+            "id": "f9b7c128-6d30-425c-87db-645ea8cda143",
            "metadata": {},
            "outputs": [],
            "source": []