{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"## Efficient zarr loading and caching\n### Brett Olsen, May 2023\n\nThis is an extension of work I did with Moshe Levy to improve the performance of Volume Annotator, but is easily made applicable to the data here.\nThe zarr loading (e.g., [here](https://www.kaggle.com/code/brettolsen/simpler-zarr-tif-image-loading)) allows easy access to slices of the full data stack but every individual data stack must be loaded separately.\nWe came up with a caching approach that allows you to hold regions *nearby* recent queries in memory, allowing extremely fast access to them.\n\nWe will still need to install the zarr package first.","metadata":{}},{"cell_type":"code","source":"!pip install zarr","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2023-05-12T16:39:01.166843Z","iopub.execute_input":"2023-05-12T16:39:01.168084Z","iopub.status.idle":"2023-05-12T16:39:17.019673Z","shell.execute_reply.started":"2023-05-12T16:39:01.168036Z","shell.execute_reply":"2023-05-12T16:39:17.018176Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import tifffile\nimport numpy as np\nimport matplotlib.pyplot as plt\nimport os\nimport time\nimport psutil\nimport zarr\nimport sys\nfrom tqdm import tqdm\nfrom collections import Counter\n\nINPUT_FOLDER = \"/kaggle/input/vesuvius-challenge-ink-detection/\"","metadata":{"execution":{"iopub.status.busy":"2023-05-12T17:07:12.043693Z","iopub.execute_input":"2023-05-12T17:07:12.044155Z","iopub.status.idle":"2023-05-12T17:07:12.051132Z","shell.execute_reply.started":"2023-05-12T17:07:12.044124Z","shell.execute_reply":"2023-05-12T17:07:12.049697Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"First I'll set up the same timing and memory tracking functions I have used previously to demonstrate the data use.","metadata":{}},{"cell_type":"code","source":"class TimerError(Exception):\n    pass\n\nclass Timer():\n    \"\"\"This is a utility class that, when used as a context manager,\n    will report the time spent on code inside its block.\n    \"\"\"\n    def __init__(self, text=None):\n        if text is not None:\n            self.text = text + \": {:0.4f} seconds\"\n        else:\n            self.text = \"Elapsed time: {:0.4f} seconds\"\n        def logfunc(x):\n            print(x)\n        self.logger = logfunc\n        self._start_time = None\n\n    def start(self):\n        if self._start_time is not None:\n            raise TimerError(\"Timer is already running.  Use .stop() to stop it.\")\n        self._start_time = time.time()\n\n    def stop(self):\n        if self._start_time is None:\n            raise TimerError(\"Timer is not running.  Use .start() to start it.\")\n        elapsed_time = time.time() - self._start_time\n        self._start_time = None\n\n        if self.logger is not None:\n            self.logger(self.text.format(elapsed_time))\n\n        return elapsed_time\n\n    def __enter__(self):\n        self.start()\n        return self\n    \n    def __exit__(self, exc_type, exc_value, exc_traceback):\n        self.stop()\n\ndef show_mem_use():\n    process = psutil.Process()\n    mb_mem = process.memory_info().rss / 1e6\n    print(f\"{mb_mem:6.2f} MB used\")","metadata":{"execution":{"iopub.status.busy":"2023-05-12T16:39:32.313951Z","iopub.execute_input":"2023-05-12T16:39:32.314368Z","iopub.status.idle":"2023-05-12T16:39:32.325124Z","shell.execute_reply.started":"2023-05-12T16:39:32.314336Z","shell.execute_reply":"2023-05-12T16:39:32.324029Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We'll go ahead and wrap the zarr loading into a single function.  Note that this function supports not just the z-indexed tiff files used in the kaggle competition but also the cubic 3D tiff cells that were reprocessed from the full scroll data.  These 3D cubes allow loading a compact region of the data by indexing into many fewer tiff files, allowing much faster access to a large number of z-slices.  If you're planning on doing a lot of work, it's well worth the time reprocessing your data into this format.\n\nIf there's sufficient demand, I'll consider adding code to reprocess the data into this format.","metadata":{}},{"cell_type":"code","source":"def load_tif(path):\n    \"\"\"This function will take a path to a folder that contains a stack of .tif\n    files and returns a concatenated 3D zarr array that will allow access to an\n    arbitrary region of the stack.\n\n    We support two different styles of .tif stacks.  The first are simply\n    numbered filenames, e.g., 00.tif, 01.tif, 02.tif, etc.  In this case, the\n    numbers are taken as the index into the zstack, and we assume that the zslices\n    are fully continuous.\n\n    The second follows @spelufo's reprocessing, and are not 2D images but 3D cells\n    of the data.  These should be labeled\n\n    cell_yxz_YINDEX_XINDEX_ZINDEX\n\n    where these provide the position in the Y, X, and Z grid of cuboids that make\n    up the image data.\n    \"\"\"\n    # Get a list of .tif files\n    tiffs = [filename for filename in os.listdir(path) if filename.endswith(\".tif\")]\n    if all([filename[:-4].isnumeric() for filename in tiffs]):\n        # This looks like a set of z-level images\n        tiffs.sort(key=lambda f: int(''.join(filter(str.isdigit, f))))\n        paths = [os.path.join(path, filename) for filename in tiffs]\n        store = tifffile.imread(paths, aszarr=True)\n    elif all([filename.startswith(\"cell_yxz_\") for filename in tiffs]):\n        # This looks like a set of cell cuboid images\n        images = tifffile.TiffSequence(os.path.join(path, \"*.tif\"), pattern=r\"cell_yxz_(\\d+)_(\\d+)_(\\d+)\")\n        store = images.aszarr(axestiled={0: 1, 1: 2, 2: 0})\n    stack_array = zarr.open(store, mode=\"r\")\n    return stack_array","metadata":{"execution":{"iopub.status.busy":"2023-05-12T16:39:34.595984Z","iopub.execute_input":"2023-05-12T16:39:34.596344Z","iopub.status.idle":"2023-05-12T16:39:34.608719Z","shell.execute_reply.started":"2023-05-12T16:39:34.596318Z","shell.execute_reply":"2023-05-12T16:39:34.604410Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Finally, we have the loader class.  I've added a couple print statements to make it more clear what's going on; feel free to comment those out\nfor production code.\n","metadata":{}},{"cell_type":"code","source":"def slice_to_hashable(slice):\n    return (slice.start, slice.stop)\n\ndef hashable_to_slice(item):\n    return slice(item[0], item[1], None)\n\nclass Loader:\n    \"\"\"Provides a cached interface to the data.\n\n    When provided with z, x, and y slices, queries whether that\n    data is available from any previously cached data.  If it\n    is, returns the subslice of that cached data that contains\n    the data we need.  \n\n    If it isn't, then we:\n    - clear out the oldest cached data if we need space\n    - load the data we need, plus some padding\n        - padding depends on the chunking of the zarr array and the\n          size of the data loaded\n    - store that data in our cache, along with the access time.\n\n    The cache is a dictionary, indexed by a namedtuple of\n    slice objects (zslice, xslice, yslice), that provides a\n    dict of \n    {\"accesstime\": last access time, \"array\": numpy array with data}\n    \"\"\"\n\n    def __init__(self, zarr_array, max_mem_gb=5):\n        self.shape = zarr_array.shape\n        print(\"Generating loader...\")\n        print(\"Data array shape: \", self.shape)\n        self.cache = {}\n        self.zarr_array = zarr_array\n        self.max_mem_gb = max_mem_gb\n        \n        full_size = self.estimate_slice_size(\n            slice(None, None, None),\n            slice(None, None, None),\n            slice(None, None, None),\n        )\n        print(f\"Estimated full array size (GB): {full_size:.2f}\")\n\n        chunk_shape = self.zarr_array.chunks\n        if chunk_shape[0] == 1:\n            self.chunk_type = \"zstack\"\n        else:\n            self.chunk_type = \"cuboid\"\n\n\n    def _check_slices(self, cache_slice, new_slice, length):\n        \"\"\"Queries whether a new slice's data is contained\n        within an older slice.\n\n        Note we don't handle strided slices.\n        \"\"\"\n        if isinstance(new_slice, int):\n            new_start = new_slice\n            new_stop = new_slice + 1\n        else:\n            new_start = 0 if new_slice.start is None else max(0, new_slice.start)\n            new_stop = length if new_slice.stop is None else min(length, new_slice.stop)\n        if isinstance(cache_slice, int):\n            cache_start = cache_slice\n            cache_stop = cache_slice + 1\n        else:\n            cache_start = 0 if cache_slice.start is None else cache_slice.start\n            cache_stop = length if cache_slice.stop is None else cache_slice.stop\n        if (new_start >= cache_start) and (new_stop <= cache_stop):\n            # New slice should be the slice in the cached data that gives the request\n            if isinstance(new_slice, int):\n                return new_start - cache_start\n            else:\n                return slice(\n                    new_start - cache_start,\n                    new_stop - cache_start,\n                    None\n                )\n        else:\n            return None\n\n    def check_cache(self, zslice, xslice, yslice):\n        \"\"\"Looks through the cache to see if any cached data\n        can provide the data we need.\n        \"\"\"\n        for key in self.cache.keys():\n            cache_zslice = hashable_to_slice(key[0])\n            cache_xslice = hashable_to_slice(key[1])\n            cache_yslice = hashable_to_slice(key[2])\n            sub_zslice = self._check_slices(cache_zslice, zslice, self.shape[0])\n            sub_xslice = self._check_slices(cache_xslice, xslice, self.shape[1])\n            sub_yslice = self._check_slices(cache_yslice, yslice, self.shape[2])\n            if sub_zslice is None or sub_xslice is None or sub_yslice is None:\n                continue\n            # At this point we have a valid slice to index into this subarray.\n            # Update the access time and return the array.\n            self.cache[key][\"accesstime\"] = time.time()\n            return self.cache[key][\"array\"][sub_zslice, sub_xslice, sub_yslice]\n        return None\n\n    @property\n    def cache_size(self):\n        total_mem_gb = 0\n        for _, data in self.cache.items():\n            total_mem_gb += data[\"array\"].nbytes / 1e9\n        return total_mem_gb\n\n    def view_cache(self):\n        for key, value in self.cache.items():\n            print(\"Access time: {}\".format(value[\"accesstime\"]))\n            print(\"Slices: {}\".format(key))\n            print(\"Est size (GB): {}\".format(value[\"array\"].nbytes / 1e9))\n    \n    def empty_cache(self):\n        \"\"\"Removes the oldest item in the cache to free up memory.\n        \"\"\"\n        if not self.cache:\n            return\n        oldest_time = None\n        oldest_key = None\n        for key, value in self.cache.items():\n            if oldest_time is None or value[\"accesstime\"] < oldest_time:\n                oldest_time = value[\"accesstime\"]\n                oldest_key = key\n        del self.cache[oldest_key]\n\n    def estimate_slice_size(self, zslice, xslice, yslice):\n        def slice_size(s, l):\n            if isinstance(s, int):\n                return 1\n            elif isinstance(s, slice):\n                s_start = 0 if s.start is None else s.start\n                s_stop = l if s.stop is None else s.stop\n                return s_stop - s_start\n            raise ValueError(\"Invalid index\")\n        return (\n            self.zarr_array.dtype.itemsize * \n            slice_size(zslice, self.shape[0]) *\n            slice_size(xslice, self.shape[1]) * \n            slice_size(yslice, self.shape[2])\n        ) / 1e9\n\n    def pad_request(self, zslice, xslice, yslice):\n        \"\"\"Takes a requested slice that is not loaded in memory\n        and pads it somewhat so that small movements around the\n        requested area can be served from memory without hitting\n        disk again.\n\n        For zstack data, prefers padding in xy, while for cuboid\n        data, prefers padding in z.\n        \"\"\"\n        def pad_slice(old_slice, length, int_add=1):\n\n            if isinstance(old_slice, int):\n                if length == 1:\n                    return old_slice\n                return slice(\n                    max(0, old_slice - int_add),\n                    min(length,old_slice + int_add + 1),\n                    None\n                )\n            start = old_slice.start if old_slice.start is not None else 0\n            stop = old_slice.stop if old_slice.stop is not None else length\n            adj_width = (stop - start) // 2 + 1\n            return slice(\n                max(0, start - adj_width),\n                min(length, stop + adj_width),\n                None\n            )\n        est_size = self.estimate_slice_size(zslice, xslice, yslice)\n        if (3 * est_size) >= self.max_mem_gb:\n            # No padding; the array's already larger than the cache.\n            return zslice, xslice, yslice\n        if self.chunk_type == \"zstack\":\n            # First pad in X and Y, then Z\n            xslice = pad_slice(xslice, self.shape[1])\n            est_size = self.estimate_slice_size(zslice, xslice, yslice)\n            if (3 * est_size) >= self.max_mem_gb:\n                return zslice, xslice, yslice\n            yslice = pad_slice(yslice, self.shape[2])\n            est_size = self.estimate_slice_size(zslice, xslice, yslice)\n            if (3 * est_size) >= self.max_mem_gb:\n                return zslice, xslice, yslice\n            zslice = pad_slice(zslice, self.shape[0], int_add=1)\n        elif self.chunk_type == \"cuboid\":\n            # First pad in Z by 5 in each direction if we have space, then in XY\n            zslice = pad_slice(\n                zslice, \n                self.shape[0], \n                int_add=min(5, self.max_mem_gb // (2 * est_size))\n            )\n            est_size = self.estimate_slice_size(zslice, xslice, yslice)\n            if (3 * est_size) >= self.max_mem_gb:\n                return zslice, xslice, yslice\n            xslice = pad_slice(xslice, self.shape[1])\n            est_size = self.estimate_slice_size(zslice, xslice, yslice)\n            if (3 * est_size) >= self.max_mem_gb:\n                return zslice, xslice, yslice\n            yslice = pad_slice(yslice, self.shape[2])\n        return zslice, xslice, yslice\n\n    def __getitem__(self, key):\n        \"\"\"Overloads the slicing operator to get data with caching\n        \"\"\"\n        zslice, xslice, yslice = key\n        print(f\"Querying {key}\")\n        for item in (zslice, xslice, yslice):\n            if isinstance(item, slice) and item.step is not None:\n                raise ValueError(\"Sorry, we don't support strided slices yet\")\n        # First check if we have the requested data already in memory\n        result = self.check_cache(zslice, xslice, yslice)\n        if result is not None:\n            print(\"Serving from cache\")\n            return result\n        # Pad out the requested slice before we pull it from disk\n        # so that we cache neighboring data in memory to avoid\n        # repeatedly hammering the disk\n        padded_zslice, padded_xslice, padded_yslice = self.pad_request(zslice, xslice, yslice)\n        print(f\"Padded query to {padded_zslice, padded_xslice, padded_yslice}\")\n        est_size = self.estimate_slice_size(padded_zslice, padded_xslice, padded_yslice)\n        print(f\"Estimated padded size (GB): {est_size:.2f}\")\n        # Clear out enough space from the cache that we can fit the new\n        # request within our memory limits.\n        while self.cache and (self.cache_size + est_size) > self.max_mem_gb:\n            self.empty_cache()\n        padding = self.zarr_array[padded_zslice, padded_xslice, padded_yslice]\n        self.cache[(\n            slice_to_hashable(padded_zslice),\n            slice_to_hashable(padded_xslice),\n            slice_to_hashable(padded_yslice),\n        )] = {\n            \"accesstime\": time.time(),\n            \"array\": padding,\n        }\n\n        result = self.check_cache(zslice, xslice, yslice)\n        if result is None:\n            # We shouldn't get cache misses!\n            print(\"Unexpected cache miss\")\n            print(zslice, xslice, yslice)\n            print(padded_zslice, padded_xslice, padded_yslice)\n            raise ValueError(\"Cache miss after cache loading\")\n        return result","metadata":{"execution":{"iopub.status.busy":"2023-05-12T16:39:35.991552Z","iopub.execute_input":"2023-05-12T16:39:35.991928Z","iopub.status.idle":"2023-05-12T16:39:36.173253Z","shell.execute_reply.started":"2023-05-12T16:39:35.991900Z","shell.execute_reply":"2023-05-12T16:39:36.171302Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Now let's try it out.  We'll load the zarr and then put that zarr array inside the loader class that handles the caching.","metadata":{}},{"cell_type":"code","source":"show_mem_use()\nwith Timer():\n    tif_path = os.path.join(INPUT_FOLDER, \"train\", \"2\", \"surface_volume\")\n    zarr_array = load_tif(tif_path)\n    data_loader = Loader(zarr_array, max_mem_gb=8)\nshow_mem_use()","metadata":{"execution":{"iopub.status.busy":"2023-05-12T16:39:39.019042Z","iopub.execute_input":"2023-05-12T16:39:39.019530Z","iopub.status.idle":"2023-05-12T16:39:41.548516Z","shell.execute_reply.started":"2023-05-12T16:39:39.019493Z","shell.execute_reply":"2023-05-12T16:39:41.547333Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"You can see that the setup was extremely fast:  a tenth of a second.  That's because we haven't actually loaded any of this data into memory yet:  the zarr array and loader class together only take a couple hundred KB of memory, even though trying to access the full array all at once would take 18GB!\n\nLet's do some queries and see how fast this is.\n\nLet's load in a whole z-slice of the fragment and take a look at it.","metadata":{}},{"cell_type":"code","source":"show_mem_use()\nwith Timer():\n    frag_slice = data_loader[35,:,:]\nplt.imshow(frag_slice)\nshow_mem_use()\ndata_loader.view_cache()","metadata":{"execution":{"iopub.status.busy":"2023-05-12T16:39:57.380215Z","iopub.execute_input":"2023-05-12T16:39:57.380618Z","iopub.status.idle":"2023-05-12T16:40:11.005983Z","shell.execute_reply.started":"2023-05-12T16:39:57.380580Z","shell.execute_reply":"2023-05-12T16:40:11.005157Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"That took a bit to load, probably mostly disk I/O.  We're using quite a bit more memory now:  the padded array we've loaded into the cache is about 850 MB, but we also have the sliced data as well.  \n\nI'm showing you the cache at the end:  the cache is indexed by access time, and when we run out of space in the cache, it will start clearing out older cached arrays to make space for new ones.\n\nLet's load up a smaller section of this same slice.","metadata":{}},{"cell_type":"code","source":"show_mem_use()\nwith Timer():\n    frag_slice = data_loader[35,5000:6000,5000:6000]\nplt.imshow(frag_slice)\nshow_mem_use()\nprint(\"Cache\")\ndata_loader.view_cache()","metadata":{"execution":{"iopub.status.busy":"2023-05-12T16:42:38.188455Z","iopub.execute_input":"2023-05-12T16:42:38.188873Z","iopub.status.idle":"2023-05-12T16:42:38.616667Z","shell.execute_reply.started":"2023-05-12T16:42:38.188842Z","shell.execute_reply":"2023-05-12T16:42:38.615703Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Now look:  our memory use went down (because we could de-allocate the much larger frag_slice we had before) and the elapsed time was ~10ms, because we could load this data directly from the cache.  \nNow if you wanted to, you could do this bit yourself because it's just a subset of the data we queried for before.\n\nHowever, because we padded our initial zarr query, we have *more* data than our initial query in the cache, which is loadable just as fast.  Let's load up something from a neighboring slice.","metadata":{}},{"cell_type":"code","source":"show_mem_use()\nwith Timer():\n    frag_slice = data_loader[34,5000:6000,5000:6000]\nplt.imshow(frag_slice)\nshow_mem_use()\nprint(\"Cache\")\ndata_loader.view_cache()","metadata":{"execution":{"iopub.status.busy":"2023-05-12T16:45:07.334862Z","iopub.execute_input":"2023-05-12T16:45:07.335310Z","iopub.status.idle":"2023-05-12T16:45:07.750048Z","shell.execute_reply.started":"2023-05-12T16:45:07.335278Z","shell.execute_reply":"2023-05-12T16:45:07.748713Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"This is still a cache hit, and still just as fast.\n\nLet's do a couple more arbitrary accesses to fill up the cache.","metadata":{}},{"cell_type":"code","source":"show_mem_use()\nwith Timer():\n    frag_slice = data_loader[30,2000:10000,2000:8000]\n    frag_slice = data_loader[25,2000:10000,2000:8000]\n    frag_slice = data_loader[42,2000:10000,2000:8000]\nplt.imshow(frag_slice)\nshow_mem_use()\nprint(\"Cache\")\ndata_loader.view_cache()","metadata":{"execution":{"iopub.status.busy":"2023-05-12T16:46:25.833194Z","iopub.execute_input":"2023-05-12T16:46:25.833641Z","iopub.status.idle":"2023-05-12T16:46:52.581191Z","shell.execute_reply.started":"2023-05-12T16:46:25.833609Z","shell.execute_reply":"2023-05-12T16:46:52.580104Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"data_loader.cache_size","metadata":{"execution":{"iopub.status.busy":"2023-05-12T16:48:39.279100Z","iopub.execute_input":"2023-05-12T16:48:39.279497Z","iopub.status.idle":"2023-05-12T16:48:39.286049Z","shell.execute_reply.started":"2023-05-12T16:48:39.279464Z","shell.execute_reply":"2023-05-12T16:48:39.285241Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"That took a bit, because we had to make a new cache access every time.  The default loader will fill up to 5GB of data, but this is adjustable on use.  Right now we're using just over 3 for the data in cache.  Let's drop this down and see what happens when we make a new query.","metadata":{}},{"cell_type":"code","source":"data_loader.max_mem_gb = 3.0","metadata":{"execution":{"iopub.status.busy":"2023-05-12T16:49:18.678997Z","iopub.execute_input":"2023-05-12T16:49:18.680012Z","iopub.status.idle":"2023-05-12T16:49:18.684373Z","shell.execute_reply.started":"2023-05-12T16:49:18.679961Z","shell.execute_reply":"2023-05-12T16:49:18.683386Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"show_mem_use()\nwith Timer():\n    frag_slice = data_loader[46,2000:10000,2000:8000]\n    frag_slice = data_loader[52,2000:10000,2000:8000]\nplt.imshow(frag_slice)\nshow_mem_use()\nprint(\"Cache\")\ndata_loader.view_cache()","metadata":{"execution":{"iopub.status.busy":"2023-05-12T16:49:35.084667Z","iopub.execute_input":"2023-05-12T16:49:35.085950Z","iopub.status.idle":"2023-05-12T16:49:58.302770Z","shell.execute_reply.started":"2023-05-12T16:49:35.085897Z","shell.execute_reply":"2023-05-12T16:49:58.301722Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"You can see that we now only have three of our slices in memory; the older ones have been dropped to make room for the new queries.\n\nLet's show how this works if we were doing some sort of systematic walk through the data.\nWe'll query a bunch of contiguous regions, storing access times for each one, and then look at the access time distributions.","metadata":{}},{"cell_type":"code","source":"# Bump up the memory use a bit\ndata_loader.max_mem_gb = 8.0\n\nclass HiddenPrints:\n    def __enter__(self):\n        self._original_stdout = sys.stdout\n        sys.stdout = open(os.devnull, 'w')\n\n    def __exit__(self, exc_type, exc_val, exc_tb):\n        sys.stdout.close()\n        sys.stdout = self._original_stdout\n\naccess_times = []\n\nfor z in range(24, 26):\n    for x in tqdm(range(3000, 3500)):\n        # We're going to hide the debug messages\n        with HiddenPrints():\n            for y in range(3000, 3500):\n                start = time.time()\n                data = data_loader[z-1:z+1, x-64:x+64, y-64:y+64]\n                stop = time.time()\n                # Yes, this is lazy and slows things down.  Won't affect the access time data, though.\n                access_times.append(stop - start)","metadata":{"execution":{"iopub.status.busy":"2023-05-12T17:00:17.565397Z","iopub.execute_input":"2023-05-12T17:00:17.565799Z","iopub.status.idle":"2023-05-12T17:06:26.136664Z","shell.execute_reply.started":"2023-05-12T17:00:17.565768Z","shell.execute_reply":"2023-05-12T17:06:26.135582Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"That took a bit:  we were doing 2 * 500 * 500 = 500,000 separate queries.  Let's look at the access time distribution.","metadata":{}},{"cell_type":"code","source":"access_times = np.array(access_times)\nplt.figure()\nplt.grid()\nplt.hist(access_times, bins=100);\nplt.yscale('log')\nplt.ylabel(\"# of accesses\")\nplt.xlabel(\"Access time (seconds)\")","metadata":{"execution":{"iopub.status.busy":"2023-05-12T17:08:43.411129Z","iopub.execute_input":"2023-05-12T17:08:43.411583Z","iopub.status.idle":"2023-05-12T17:08:44.126259Z","shell.execute_reply.started":"2023-05-12T17:08:43.411548Z","shell.execute_reply":"2023-05-12T17:08:44.125039Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"You can see the distribution is bimodal: most accesses are very fast (because we just pulled from the cache) while a small percentage had to do a fresh query and go to disk to get the data.  Let's quantify that a bit.","metadata":{}},{"cell_type":"code","source":"threshold = 0.4\nfast = access_times < threshold\nslow = access_times >= threshold\nprint(fast.sum(), slow.sum())\nprint(access_times[fast].mean(), access_times[slow].mean())","metadata":{"execution":{"iopub.status.busy":"2023-05-12T17:10:48.708669Z","iopub.execute_input":"2023-05-12T17:10:48.709191Z","iopub.status.idle":"2023-05-12T17:10:48.722521Z","shell.execute_reply.started":"2023-05-12T17:10:48.709152Z","shell.execute_reply":"2023-05-12T17:10:48.721182Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"64 cache misses across all half a million queries, which took ~0.7 seconds each, while the cache hits were a thousand times faster, at ~0.6ms.\n\nNow, this is obviously overkill for this particular setup:  you could just hold all the data you're working with here in memory at once.  But for the larger fragments, and most especially for the full scroll data (which is ~8TB), it is extremely useful to have easy and fast access to contiguous regions of the data.","metadata":{}}]}