{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"toc": true
},
"source": [
"Requires **HyperSpy 1.4.2 or above**\n",
"\n",
"This tutorial introduces to the processing of large dataset - which can not fit into memory - using HyperSpy. It introduce the concept of out-of-core computation algorithms (also refer as lazy processing) and the main difference between lazy and non-lazy processing as well as technicallities you need to be aware of to optimise performance.\n",
"The corresponding section of the HyperSpy documentation is [the big data section](https://hyperspy.readthedocs.io/en/stable/user_guide/big_data.html#limitations).\n",
"\n",
"### Credits and changes\n",
"\n",
"* 29/07/2019 Eric Prestat. Add more details and introduction for the M&M Sunday short course.\n",
"* 15/03/2019 Francisco de la Peña. Create tutorial for the HyperSpy workshop at ePSIC.\n",
"\n",
"## Table of contents\n",
"1. [Introduction to lazy processing](#1.-Introduction-to-lazy-processing)\n",
"2. [Loading data lazily](#2.-Loading-data-lazily)\n",
"3. [Plotting lazily](#3.-Plotting-lazily)\n",
"4. [Rebinning](#4.-Rebinning)\n",
"5. [ROI in navigation dimension](#5.-ROI-in-navigation-dimension)\n",
"6. [Summary](#6.-Summary)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 1. Introduction to lazy processing\n",
"\n",
"Lazy processing refers to the use of [out-of-core computation algorithms](https://en.wikipedia.org/wiki/External_memory_algorithm) to process very large data, which are usually too large to fit into the computer's memory at one time. The main idea is to chunk the data in pieces, small enough, that can be processed in memory as illustrated by the following diagram:\n",
"\n",
"\n",
"\n",
"HyperSpy internally uses the [dask library](https://docs.dask.org/en/latest/index.html), which implements the numpy interface to larger-than-memory or distributed environments. The typically workflow for processing data lazily is available on a disk:\n",
" 1. \"load\" data from disk with a defined chunking\n",
" 2. schedule operations\n",
" 3. do the computation\n",
"\n",
"**Steps 1 and 2 are very fast**, because nothing is actually done, other than initialising and scheduling the tasks to be peformed.\n",
"**Step 3 is slow**, because all the computation is performed at this stage. Most of the time, this is signficantly slower than in memory processing, because the chunks of data needs to be read and written from/to disk when on request of the scheduler.\n",
"\n",
"The following diagram shows a task graph, where square and rounds represent arrays and functions, respectively. This graph is an example of how dask is going process the data from the large data set on the harddrive (HDD) into chunks and do computation on these.\n",
"\n",
"This example below a simple example shows of how to perform out-of-core computation using dask, here the calculation of the sum of an array:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import dask.array as da\n",
"# Create a 15x15 array filled with ones and chunks size to 5x5\n",
"x = da.ones((15, 15), chunks=(5, 5))\n",
"\n",
"# Take the sum()\n",
"y = x.sum()\n",
"\n",
"# do the computation\n",
"y.compute()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The corresponding tasks can be representated by the following graph (square are array):\n",
" - square are arrays, which can fit into memory\n",
" - circle are operations\n",
" \n",
"\n",
"\n",
"Read [the graph](https://docs.dask.org/en/latest/graphs.html) section of the dask documenation for more explanation. The taks graph can be visualised using the [visualise](https://docs.dask.org/en/latest/graphviz.html) method.\n",
"\n",
"For a short ontroductory presentation on dask and its principle see http://matthewrocklin.com/slides/plotcon-2016.html. The following animation - taken from the previously mentioned presentation - illustrates the execution of the scheduled computation of many tasks:\n",
"\n",
"![Dask's directed acyclic diagram (DAG)](lazy/grid_search_schedule.gif)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The implementation of out-of-core computation in HyperSpy aims to make processing very large data (not fitting into memory) as seamlessly as possible and similar to in-memory data. This tutorial covers the main difference between lazy and non-lazy processing as well as technicallities you need to be aware of to optimise performance.\n",
"\n",
"## 2. Loading data lazily\n",
"\n",
"As usual, we start by setting up the matplotlib backend and importing hyperspy"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib qt\n",
"import hyperspy.api as hs"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's start by setting the ``logging level`` to \"INFO\" using the preferences GUI as below. Once done, click ``Save`` and restart the kernel."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Open the preferences GUI and set the logging level to INFO\n",
"hs.preferences.gui()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"For this tutorial we are going to start by loading a large spectrum image lazily that won't fit in the RAM of a standard laptop computer in 2018, so notice the ``lazy`` keyword and don't forget to set it to ``True``."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Open lazily the file \"lazy_demo_data_hyperspy_chunking.hspy\" in the subfolder \"lazy\":\n",
"s = hs.load(\"lazy/lazy_demo_data_hyperspy_chunking.hspy\", lazy=True)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's check what sort of object we have stored in the ``s`` variable"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"print(s)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This is a spectrum image with *4096* spectral channels and *512x512* image dimensions. Its size in GB:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Use the \"nbytes\" attribute of the numpy array to calculate the size on disk\n",
"print(s.data.nbytes / 1e9)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"That's more than *8 GB* of data. However, loading it tooks no time. That's because HyperSpy didn't load the data yet. It'll do it on demand.\n",
"\n",
"## 3. Plotting lazily \n",
"\n",
"Let's plot the spectrum image. Usually we would call simply:\n",
"\n",
"```python\n",
"s.plot()\n",
"```\n",
"\n",
"And a navigator would appear. However, in this case, computing the navigator would take too long. Therefore, we start by navigating the dataset \"blindly\" using sliders.\n",
"\n",
"Notice the ``Continous update`` checkbox in the ipywidgets (in you have the ``hyperspy_gui_ipywidgets`` installed and well configured). Unchecking it won't update the signal figure until you release the slider, what is handy to avoid lagging when dealing with large datasets."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Plot the `s` signal by setting \"slider\" to the parameter `navigator`:\n",
"s.plot(navigator=\"slider\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"If we were not working in lazy mode we could easily explore this dataset as an image stack by transposing and plotting it as follows:\n",
"\n",
"```python\n",
"s.T.plot(navigator=\"slider\")\n",
"```\n",
"\n",
"However, executing the code above with the current file risks saturating your computer memory. This is because of the way the data is chunked in the file.\n",
"\n",
"Let's have a look at the data chunks:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Get the chunks information by calling the `print` function of `s.data`:\n",
"print(s.data)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This means that the data is stored in the file as *4x4x4096* blocks with the biggest (4096) dimension corresponding to the spectral dimension. In other words, the chunks are *4 px x 4 px* spectrum images. This chunking is good if we want to perform operations that operate on indiviual spectra. However, it is a terrible chunking configuration for e.g. slicing the spectrum image as this operation would require loading the whole dataset from the disk.\n",
"\n",
"Let's load the same dataset but stored with a different chunking, in this case optimised for slicing it:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Open lazily the file \"lazy_demo_data_hyperspy_image_chunking.hspy\" in the subfolder \"lazy\":\n",
"im = hs.load(\"lazy/lazy_demo_data_hyperspy_image_chunking.hspy\", lazy=True)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's check what sort of object we've loaded"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"print(im)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"It's the same dataset, but now the \"view\" is different: it is configured as an image stack.\n",
"\n",
"We could have obtained the same view by transposing the spectrum image as follow"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"sT = s.T\n",
"print(sT)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"It looks the same, but when we look at if we look at the chunks we'll see that they're very different:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"print(sT.data)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"For the transposed dataset the chunks are just like in the file, good for looking at spectra but bad for looking at images.\n",
"\n",
"Let's have a look at the chunking of the ``lazy_demo_data_hyperspy_image_chunking.hspy`` file:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"print(im.data)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In this case, the chunks are individual images, what's optimal for slicing the spectrum image."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"im.plot(navigator=\"slider\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We notice that navigation the image stack is fast, despite the fact that the data is not all loaded into memory at once but only 1 image at a time on demand.\n",
"\n",
"### 3.1 Creating a \"navigator\" for the two views of the dataset\n",
"\n",
"Let's create a \"navigator\" for the spectrum image the way that HyperSpy usually do it, but summing all the signal dimensions.\n",
"\n",
"This operation could be performed as follows for the spectrum image:\n",
"\n",
"```python\n",
"# by default, the sum is performed over navigation axes, no need to specify the axis argument here.\n",
"s_navigator = s.sum()\n",
"```\n",
"\n",
"However, in this case, due to the particular chunking of this file, we know that this is a bad idea as it'll require loading the whole dataset in order to perform the computation. Instead we could sum all the navigation dimensions of the image stack which is equivalent but a lot more efficient in this case."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Calculate the sum over the navigation axes using the `sum` method of `im`:\n",
"s_navigator = im.sum()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Notice the log messages informing us that HyperSpy has automatically rechunked the dataset. This is in order to optimize the chunks for the operation requested. However, this rechunking is \"virtual\" i.e. HyperSpy obviously didn't modify the chunks in the file. This is still useful to define the size of the chunks that will be processed by the ``sum`` function.\n",
"\n",
"Notice also that performing the sum took no time. That's because so far it hasn't actually performed the operation. We can confirm this by verifying that the object is a ``LazySignal``"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"print(s_navigator)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In order to compute the sum we have to call the ``compute`` method."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Create the s_navigator by summing over the navigation axes:\n",
"s_navigator = im.sum()\n",
"# Compute the lazy signal `s_navigator`\n",
"s_navigator.compute()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"print(s_navigator)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now ``s_navigator`` is no longer a ``LazySignal``.\n",
"\n",
"We can generate a navigator for the image stack in the same way:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Do the same for `im_navigator`:\n",
"im_navigator = s.sum()\n",
"# Compute the lazy signal `im_navigator`\n",
"im_navigator.compute()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's now use the navigators to navigate the spectrum image and image stack"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# plot the `s` lazy signal using the `s_navigator` as navigator\n",
"s.plot(navigator=s_navigator)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# plot the `im` lazy signal using the `im_navigator` as navigator\n",
"im.plot(navigator=im_navigator)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 4. Rebinning"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Load the data \"lazy_demo_data_hyperspy_image_chunking.hspy\" from the subfolder \"lazy\"\n",
"im = hs.load(\"lazy/lazy_demo_data_hyperspy_image_chunking.hspy\", lazy=True)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Get information about the shape and the chunks of the dara using the `print` function:\n",
"print(im.data)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We have 4096 chunks of size (1, 512, 512)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Rechunk the data automatically \n",
"im._make_lazy(rechunk=\"dask_auto\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Check again the information of the array to compare with the previous chunking\n",
"print(im.data)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now, we have significantly less chunks: 64 (4096/64) chunks instead of 4096 and there are all of egal size (64, 512, 512)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Rebin the signal dimention to 128x128 using the `rebin` method (you may need to check the documentation of `rebin`)\n",
"srb = im.rebin((4096, 128, 128))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Call `compute` to perform the computation\n",
"srb.compute()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The results of the `rebin` function is not a lazy signal"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Plot the result, which have been computed and is now in memory\n",
"srb.plot()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 5. ROI in navigation dimension"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Plot `s` with precomputed navigator image `s_navigator`\n",
"s.plot(navigator=s_navigator)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's create a rectagular ROI"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Create a RectangularROI defined by the position (left=100, top=100, right=110, bottom=110)\n",
"s_roi = hs.roi.RectangularROI(left=100, top=100, right=110, bottom=110)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Check the ROI\n",
"print(s_roi)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Calling this object extracts a rectangular roi from a HyperSpy signal. For example:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Call the ROI from a signal:\n",
"print(s_roi(s))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can change its parameters interactively displaying the associated widget as follows:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Call the `gui` method of the ROI to display the widget:\n",
"s_roi.gui()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Check the changes in the roi:\n",
"print(s_roi(s))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"And also by adding a widget to the navigator plot"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Add the ROI to the signal using the `add_widget` method of the roi:\n",
"w = s_roi.add_widget(s)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"So far, we have to execute manually the ROI on the spectrum image every time that we change its parameters. We can perform this operation automatically instead:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Use the `hs.interactive` function to automatically calculate a new signal from `s` and defined by the ROI:\n",
"s2 = hs.interactive(s_roi, signal=s,\n",
" event=s_roi.events.changed,\n",
" recompute_out_event=None)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Get the information about `s2`:\n",
"print(s2)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can use the ``interactive`` function to compute the sum of the extracted spectrum image interactively as follows:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Use the `hs.interactive` function to automatically calculate the sum of `s2`:\n",
"s2_sum = hs.interactive(s2.sum, rechunk=False,\n",
" event=s2.events.data_changed,\n",
" recompute_out_event=None)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Check the information of `s2_sum`:\n",
"s2_sum"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"And plot it. Notice that the spectrum should update when we change the ROI parameters"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"s2_sum.plot()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"What if we want to integrate the peaks in the spectra and display the corresponding image?\n",
"\n",
"When working with in-memory data we could do it using the spectrum image but, in lazy mode, it is better to work with the image stack version of the dataset."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"im.plot(navigator=im_navigator)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"im_roi = hs.roi.SpanROI(100, 200)\n",
"im_w = im_roi.add_widget(im)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"im_roi.gui()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"im2 = hs.interactive(im_roi, signal=im,\n",
" event=im_roi.events.changed,\n",
" recompute_out_event=None)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"im2_sum = hs.interactive(im2.sum, rechunk=False,\n",
" event=im2.events.data_changed,\n",
" recompute_out_event=None)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"im2_sum.plot()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# 6. Summary\n",
"\n",
"Most operations can be performed *lazily* in HyperSpy:\n",
"1. Visualisation\n",
"2. Slicing and indexing\n",
"3. Generic mathematical operations\n",
"4. Machine learning\n",
"5. Curve fitting\n",
"\n",
"See [the big data section](https://hyperspy.readthedocs.io/en/stable/user_guide/big_data.html#limitations) of the HyperSpy documentation for more information and to learn about the main difference between lazy and non-lazy signal."
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.6"
},
"toc": {
"nav_menu": {},
"number_sections": true,
"sideBar": true,
"skip_h1_title": false,
"title_cell": "Table of Contents",
"title_sidebar": "Contents",
"toc_cell": true,
"toc_position": {},
"toc_section_display": true,
"toc_window_display": false
}
},
"nbformat": 4,
"nbformat_minor": 4
}