This tutorial shows the basics of model (currently only 1D) fitting in HyperSpy from the grounds up. All data is artificial and can be generated by running the code at the end of the notebook.
13/04/2015 Tomas Ostasevicius - Developed for HyperSpy workshop at University of Cambridge
01/06/2016 Tomas Ostasevicius - updated and expanded for HyperSpy workshop at Scandem conference 2016
22/07/2016 Tomas Ostasevicius - updated for HyperSpy version 1.0
In order to use fitting in HyperSpy more effectively, it is useful to understand our structure for curve fitting.
There are three main things, related to fitting:
1. Model can be thought of as a simple box (cooking pot), where we have to put our ingredients. Without anything inside, it is not of much use in this case. Once we add some things to it and mix it a bit (do the actual fitting), however, we have our complete dish!
2. Component is the main building block (ingredient) of our model. Here we mix and match what components we need (or want) for the particular case of signal.
Each of the components is ultimately just a function that has variables that change the (shape of the) output. Such a variable in HyperSpy is called a parameter
3. Parameter is the knob that the fitting routine adjusts for a good fit. Each component must include at least one parameter in order to be able to change when fitting. A parameter is also the object that we may limit or have to adjust when the result of the fit is not satisfactory.
Ultimately, a parameter is the only important thing, as far as the fitting is concerned - components are just smart and convenient boxes to combine parameters into functions, and a model is just a box for a collection of components.
For now, let's just keep the rough structure in our heads and look at other things!
%matplotlib nbagg import hyperspy.api as hs
WARNING:hyperspy_gui_traitsui:The nbAgg matplotlib backend is not supported by the installed traitsui version and the ETS toolkit has been set to null. To set the ETS toolkit independently from the matplotlib backend, set it before importing matplotlib. WARNING:hyperspy_gui_traitsui:The traitsui GUI elements are not available.
Once imported, all the HyperSpy commands are available via the
interface. You can also look for the help with any Python object like this
Help on module hyperspy.api in hyperspy: NAME hyperspy.api - All public packages, functions and classes are available in this module. DESCRIPTION When starting HyperSpy using the ``hyperspy`` script (e.g. by executing ``hyperspy`` in a console, using the context menu entries or using the links in the ``Start Menu``, the :mod:`~hyperspy.api` package is imported in the user namespace as ``hs``, i.e. by executing the following: >>> import hyperspy.api as hs (Note that code snippets are indicated by three greater-than signs) We recommend to import the HyperSpy API as above also when doing it manually. The docstring examples assume that `hyperspy` has been imported as `hs`, numpy as ``np`` and ``matplotlib.pyplot`` as ``plt``. Functions: create_model Create a model for curve fitting. get_configuration_directory_path Return the configuration directory path. load Load data into BaseSignal instances from supported files. preferences Preferences class instance to configure the default value of different parameters. It has a CLI and a GUI that can be started by execting its `gui` method i.e. `preferences.gui()`. stack Stack several signals. interactive Define operations that are automatically recomputed on event changes. set_log_level Convenience function to set HyperSpy's the log level. The :mod:`~hyperspy.api` package contains the following submodules/packages: :mod:`~hyperspy.api.signals` `Signal` classes which are the core of HyperSpy. Use this modules to create `Signal` instances manually from numpy arrays. Note that to load data from supported file formats is more convenient to use the `load` function. :mod:`~hyperspy.api.model` Contains the :mod:`~hyperspy.api.model.components` module with components that can be used to create a model for curve fitting. :mod:`~hyperspy.api.eds` Functions for energy dispersive X-rays data analysis. :mod:`~hyperspy.api.material` Useful functions for materials properties and elements database that includes physical properties and X-rays and EELS energies. :mod:`~hyperspy.api.plot` Plotting functions that operate on multiple signals. :mod:`~hyperspy.api.datasets` Example datasets. :mod:`~hyperspy.api.roi` Region of interests (ROIs) that operate on `BaseSignal` instances and include widgets for interactive operation. :mod:`~hyperspy.api.samfire` SAMFire utilities (strategies, Pool, fit convergence tests) For more details see their doctrings. DATA preferences = <hyperspy.defaults_parser.Preferences object> FILE /home/fjd29/Python/hyperspy3/hyperspy/api.py
Help on module hyperspy.components1d in hyperspy: NAME hyperspy.components1d - Components that can be used to define a 1D model for e.g. curve fitting. DESCRIPTION There are some components that are only useful for one particular kind of signal and therefore their name are preceded by the signal name: eg. eels_cl_edge. Writing a new template is really easy, just edit _template.py and maybe take a look to the other components. For more details see each component docstring. ==================================================================== Arctan Arctan function component.. Bleasdale Bleasdale function component... DoublePowerLaw .. EELSCLEdge EELS core loss ionisation edge from hydroge.. Erf Error function component.. Exponential Exponentian function components.. Expression Create a component from a string expression.. Gaussian Normalized gaussian function component.. GaussianHF Normalized gaussian function component, wit.. HeavisideStep The Heaviside step function.. Logistic Logistic function component.. Lorentzian Cauchy-Lorentz distribution (a.k.a. Lorentz.. Offset Component to add a constant value in the y-.. PESCoreLineShape .. Polynomial n-order polynomial component... PowerLaw Power law component.. RC .. SEE Secondary electron emission component for P.. ScalableFixedPattern Fixed pattern component with interpolation .. Vignetting .. Voigt Voigt profile component with support for sh.. VolumePlasmonDrude Drude volume plasmon energy loss function c.. FILE /home/fjd29/Python/hyperspy3/hyperspy/components1d.py
First you should have a spectrum (a particular kind of the
Signal subclass!) you want to fit. Let's load a synthetic dataset with some curves named
and have a look at it.
If you can't load the dataset, it means you most likely have not generated it yet. Please run the two cells at the end of the notebook to do so.
s = hs.load("two_peaks.hspy")
<Signal1D, title: Two gaussians, dimensions: (|1024)>
Creating a model now is simple - just pass the spectrum to the function
model_reference = signal_reference.create_model()
Let's reference the model by "
m = s.create_model()
Let's look what's inside:
# | Attribute Name | Component Name | Component Type ---- | ------------------- | ------------------- | -------------------
As we can see, the model is still empty. That will not always be the case - for some types of signals, an automatic background component is added when creating a model, hence it's always good to check.
We can plot the model in exactly the same way as the signal:
The only difference from the model plot is that each data point is displayed individually.
To do anything with the model, we should create some components and add them. Let's create two gaussians, referenced as "g1" and "g2":
P.S.: keep in mind that creating a component is a function - hence there should be brackets at the end! Such as
our_component_reference = hs.model.components1D.example_component()
g1 = hs.model.components1D.Gaussian() g2 = hs.model.components1D.Gaussian()
... and add the components to our model. For that there are generally two ways:
or in lists (i.e. grouped by square brackets)
Let's check how the model looks now:
# | Attribute Name | Component Name | Component Type ---- | ------------------- | ------------------- | ------------------- 0 | Gaussian | Gaussian | Gaussian 1 | Gaussian_0 | Gaussian_0 | Gaussian
For our convenience we can rename the components as we choose, for example "large" and "small" (note that the "g1" and "g2" are only references we created for them, not names of the components)
g1.name = "large" g2.name = "small"
We can look at the model again to see the result
# | Attribute Name | Component Name | Component Type ---- | ------------------- | ------------------- | ------------------- 0 | large | large | Gaussian 1 | small | small | Gaussian
Components Parameter Value large A 1 centre 0 sigma 1 small A 1 centre 0 sigma 1
To access the values, we have to look inside the components for the parameters. It can simply be done by following the pattern:
In this case the component references are the g1 and g2, while parameter names are centre, A and sigma.
g1.sigma.value = 30
Components Parameter Value large A 1 centre 0 sigma 30 small A 1 centre 0 sigma 1
For convenience, we can also set values "in bulk" for all components in the model. The required command is
Set the area ("A" parameter) of both peaks to 500
Components Parameter Value large A 500 centre 0 sigma 30 small A 500 centre 0 sigma 1
When using HyperSpy in a notebook (like this one), we can also use additional functions to adjust parameters with a mouse (interactively).
In order to enable it, call
.notebook_interaction() for either model, component or parameter directly:
Finally, let's fit the our model and plot it afterwards to see how well (or poorly) we did
Components Parameter Value large A 29818.6 centre 52.3156 sigma 29.9525 small A 535.558 centre 71.2146 sigma 3.21026
s2 = hs.load('smoothly_moving_peaks.hspy') m2 = s2.create_model()
We can still use the two components we had in the previous part (
g2), so let's just add them to the new model and plot it to see what it looks like:
We can see that the small peak is not where it should be, so just adjust its position as previously
Now we can perform a fit for all of the signals sequentially by passing command
Even though it slows HyperSpy down, it is useful (for now!) to have the model plot windows open when the command is passed:
We can see the total result of the fit, however it is often useful to see each component individually
Now that we have the artificial spectral image (SI) fitted, we can look at it however we want! For example, when fitting, you might have noticed that the position of the small peak shifted from pixel to pixel. We can plot the centre parameter of the component to have a look at it:
Once the fit was performed, chi-squared ($\chi^2$), degrees of freedom and reduced chi-squared of the fit are automatically calculated.
They are accessible with, respectively:
Let's have a look at reduced $\chi^2$ by plotting it
We can see that the map looks relatively uniform (albeit with some noise, as expected). The colorbar scale tells us that all of the points fall in the (0.9, 1.15) range, which shows a good fit.
We can of course plot the histogram of the same data like so:
Lets say we have a slightly stranger signal that we want to fit, like this one:
s = hs.load('wobbly_peak.hspy')
It's (as the name implies) composed of a sinus + gaussian + 2nd degree polynomial. However we don't have a
sin component in the in-build library, so we'll just write our own:
sin = hs.model.components1D.Expression('A*sin(b*x + c)', name='sin',)
Then just create and add all the additional components we might need: a gaussian and a polynomial
m = s.create_model() gaus = hs.model.components1D.Gaussian() poly = hs.model.components1D.Polynomial(2) m.extend([sin, gaus, poly])
Components Parameter Value sin A 0 b 0 c 0 Gaussian A 1 centre 0 sigma 1 Polynomial coefficients 0 coefficients 0 coefficients 0
The initial values do not seem to be very useful, so let's just plot the model, turn on the widgets, and we'll play until things seem close enough:
And then fit it and look at the results!
Components Parameter Value sin A -3 b 0.100033 c -0.707963 Gaussian A 2999 centre 149.95 sigma 49.9833 Polynomial coefficients -2.19049e-11 coefficients 0.06002 coefficients 3
from scipy.ndimage import gaussian_filter import numpy as np import hyperspy.api as hs # set the parameters: blurs = [0., 1.] radius = 3 #radius of a different region domain = 10 #size of the square domain small_centres = (500, 900) small_amplitudes = (8000, 5000) # work: total = None cent = (domain//2, domain//2) y,x = np.ogrid[-cent:domain-cent, -cent:domain-cent] mask = x*x + y*y <= radius*radius for blur in blurs: s = hs.signals.Signal1D(np.ones((domain,domain, 1024))) cent = tuple([int(0.5*i) for i in s.data.shape[:-1]]) m0 = s.create_model() gs01 = hs.model.components1D.Gaussian() m0.append(gs01) gs01.sigma.value = 300 gs01.centre.map['values'][:] = (np.random.random((domain,domain)) - 0.5)*50 + small_centres gs01.centre.map['is_set'][:] = True gs01.A.map['values'][:] = 1000 * np.random.random((domain,domain)) + 300000 gs01.A.map['is_set'][:] = True gs02 = hs.model.components1D.Gaussian() m0.append(gs02) gs02.sigma.value = 30 gs02.centre.map['values'][:] = (np.random.random((domain,domain)) - 0.5)*50 + small_centres gs02.centre.map['values'][mask] = (np.random.random(gs02.centre.map['values'][mask].shape) - 0.5)*50 + 700 #gs02.centre.map['values'][10:20,10:20] = (np.random.random((10,10)) - 0.5)*100 + 200 gs02.centre.map['values'] = gaussian_filter(gs02.centre.map['values'], blur) gs02.centre.map['is_set'][:] = True gs02.A.map['values'][:] = small_amplitudes gs02.A.map['values'][mask] = small_amplitudes gs02.A.map['values'] = gaussian_filter(gs02.A.map['values'], blur) gs02.A.map['is_set'][:] = True s11 = m0.as_signal() if total is None: total = s11.data.copy() else: total = np.concatenate((total, s11.data), axis=1) s = hs.signals.Signal1D(total) s.metadata.General.author = 'Tomas Ostasevicius' s.metadata.General.title = 'Two gaussians' s.add_poissonian_noise() s.estimate_poissonian_noise_variance() s.axes_manager.name = "x" s.axes_manager.units = "nm" s.axes_manager.name = "y" s.axes_manager.units = "nm" s.axes_manager.name = "Energy" s.axes_manager.name = "eV" s.axes_manager.scale = 0.1 s.inav[5,5].save("two_peaks", overwrite=True) s.inav[10:20,:].save("smoothly_moving_peaks", overwrite=True)
import numpy as np import hyperspy.api as hs k = 1 alpha = 15 amp = 3 gaus_position = 15 gaus_width = 5 gaus_A = 300 gradient = 0.6 offset= 3 sin_component = hs.model.components1D.Expression('A * sin(k*x + alpha)', name='sin', k=k, alpha=alpha, A=amp) gaus = hs.model.components1D.Gaussian(A=gaus_A, sigma=gaus_width, centre=gaus_position) poly = hs.model.components1D.Polynomial(1) poly.coefficients.value = (gradient, offset) axis = np.linspace(0, 30, 3000, dtype='double') result = sin_component.function(axis)+ gaus.function(axis) + poly.function(axis) s = hs.signals.Signal1D(result) s.axes_manager.name = 'x' s.axes_manager.scale = 0.1 s.axes_manager.offset = 0 s.metadata.General.author = 'Tomas Ostasevicius' s.metadata.General.title = 'Sin + poly(2) + Gaussian' s.save('wobbly_peak', overwrite=True)