HyperSpy Fitting tutorial

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.

Required version:

HyperSpy 1.3

Author:

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

Terminology and relationships

[back to top]

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.

Examples:

  • Lorentzian (Cauchy)
  • Gaussian
  • Voigt (a combination of Lorentzian and Gaussian)
  • Offset (i.e. constant background)
  • Exponential function
  • ...
  • [create your own or use the very specialised ones!]

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!

Importing

[back to top]

HyperSpy, like many other Python libraries, first has to be imported in your Python setup in order to be used. Once it is, all the relevant commands can be looped up using the autocompletion feature of the IPython.

Lets import the HyperSpy and set up plotting.

In [1]:
%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.

How do I .... again?

[back to top]

Once imported, all the HyperSpy commands are available via the

hs.<something>

interface. You can also look for the help with any Python object like this

help(<something>)

In [2]:
help(hs)
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


In [3]:
help(hs.model.components1D)
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


Fitting in HyperSpy (detailed)

[back to top]

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

"two_peaks.hspy"

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.

In [4]:
s = hs.load("two_peaks.hspy")
In [5]:
s
Out[5]:
<Signal1D, title: Two gaussians, dimensions: (|1024)>
In [6]:
s.plot()

[back to top]

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".

In [7]:
m = s.create_model()

Let's look what's inside:

In [8]:
m.components
Out[8]:
   # |      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:

In [9]:
m.plot()

The only difference from the model plot is that each data point is displayed individually.

[back to top]

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()

In [10]:
g1 = hs.model.components1D.Gaussian()
g2 = hs.model.components1D.Gaussian()

[back to top]

... and add the components to our model. For that there are generally two ways:

Individually

our_model_reference.append(our_component_reference)

or in lists (i.e. grouped by square brackets)

our_model_reference.extend([first_component_reference, second_component_reference])

In [11]:
m.extend([g1, g2])

Let's check how the model looks now:

In [12]:
m.components
Out[12]:
   # |      Attribute Name |      Component Name |      Component Type
---- | ------------------- | ------------------- | -------------------
   0 |            Gaussian |            Gaussian |            Gaussian
   1 |          Gaussian_0 |          Gaussian_0 |            Gaussian

[back to top]

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)

In [13]:
g1.name = "large"
g2.name = "small"

We can look at the model again to see the result

In [14]:
m.components
Out[14]:
   # |      Attribute Name |      Component Name |      Component Type
---- | ------------------- | ------------------- | -------------------
   0 |               large |               large |            Gaussian
   1 |               small |               small |            Gaussian

[back to top]

To finally see the full structure (the one we looked at here), we can print all of the parameter values of all components of the model.

In [15]:
m.print_current_values()
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:

some_component_reference.parameter_name.value

In this case the component references are the g1 and g2, while parameter names are centre, A and sigma.

In [16]:
g1.sigma.value
Out[16]:
1.0

[back to top]

We can set parameter values in exactly the same way. Let set g1 sigma value to 30:

In [17]:
g1.sigma.value = 30
In [18]:
m.print_current_values()
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

m.set_parameters_value

Set the area ("A" parameter) of both peaks to 500

In [19]:
m.set_parameters_value('A', 500)
In [20]:
m.print_current_values()
Components	Parameter	Value
large
		A	500
		centre	0
		sigma	30
small
		A	500
		centre	0
		sigma	1

[back to top]

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:

In [21]:
m.plot()
In [22]:
m.gui()

[back to top]

Finally, let's fit the our model and plot it afterwards to see how well (or poorly) we did

In [23]:
m.fit()
In [24]:
m.print_current_values()
Components	Parameter	Value
large
		A	29818.6
		centre	52.3156
		sigma	29.9525
small
		A	535.558
		centre	71.2146
		sigma	3.21026

Fitting many signals

[back to top]

When the signal is a collection of spectra (e.g. EELS / EDS), we can fit all of them automatically (sequentially).

First load the spectrum smoothly_moving_peaks.hspy, reference it as s2 and create a new model for it m2:

In [25]:
s2 = hs.load('smoothly_moving_peaks.hspy')
m2 = s2.create_model()

We can still use the two components we had in the previous part ( g1 and g2), so let's just add them to the new model and plot it to see what it looks like:

In [26]:
m2.extend([g1,g2])
m2.plot()

We can see that the small peak is not where it should be, so just adjust its position as previously

In [27]:
m2.gui()

Now we can perform a fit for all of the signals sequentially by passing command

m2.multifit()

Even though it slows HyperSpy down, it is useful (for now!) to have the model plot windows open when the command is passed:

  • In the signal widow you can see (however briefly) all of the signals and the corresponding fit
  • In the navigator window you can track the order the signals are being fitted. If you end up using this feature often, you'll notice that it has a huge impact on the accuracy of the results if chosen poorly. The next HyperSpy version (1.0) will have an automatic routine "SAMFire" to do that in a smart way, hence enabling much more complex models to be fitted effortlessly.
In [28]:
m2.multifit()

[back to top]

We can see the total result of the fit, however it is often useful to see each component individually

In [29]:
m2.plot(plot_components=True)

[back to top]

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:

In [30]:
g2.centre.plot()

Statistical analysis of fitting results

[back to top]

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:

m2.chisq

m2.dof

m2.red_chisq

Let's have a look at reduced $\chi^2$ by plotting it

In [31]:
m2.red_chisq.plot()

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:

In [32]:
m2.red_chisq.get_histogram().plot()

More interesting components

[back to top]

Lets say we have a slightly stranger signal that we want to fit, like this one:

In [33]:
s = hs.load('wobbly_peak.hspy')
In [34]:
s.plot()

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:

In [35]:
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

In [36]:
m = s.create_model()
gaus = hs.model.components1D.Gaussian()
poly = hs.model.components1D.Polynomial(2)
m.extend([sin, gaus, poly])
In [37]:
m.print_current_values()
Components	Parameter	Value
sin
		A	0
		b	0
		c	0
Gaussian
		A	1
		centre	0
		sigma	1
Polynomial
		coefficients[0]	0
		coefficients[1]	0
		coefficients[2]	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:

In [38]:
m.plot()

m.gui()

And then fit it and look at the results!

In [39]:
m.fit()
In [40]:
m.print_current_values()
Components	Parameter	Value
sin
		A	-3
		b	0.100033
		c	-0.707963
Gaussian
		A	2999
		centre	149.95
		sigma	49.9833
Polynomial
		coefficients[0]	-2.19049e-11
		coefficients[1]	0.06002
		coefficients[2]	3

Generating synthetic data

Two peaks

[back to top]

In [ ]:
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[0]:domain-cent[0], -cent[1]:domain-cent[1]]
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[0]
    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[1]

    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[0]
    gs02.A.map['values'][mask] = small_amplitudes[1]
    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[0].name = "x"
s.axes_manager[0].units = "nm"
s.axes_manager[1].name = "y"
s.axes_manager[1].units = "nm"

s.axes_manager[2].name = "Energy"
s.axes_manager[2].name = "eV"
s.axes_manager[2].scale = 0.1

s.inav[5,5].save("two_peaks", overwrite=True)
s.inav[10:20,:].save("smoothly_moving_peaks", overwrite=True)

Wobbly peak

[back to top]

In [7]:
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[0].name = 'x'
s.axes_manager[0].scale = 0.1
s.axes_manager[0].offset = 0

s.metadata.General.author = 'Tomas Ostasevicius'
s.metadata.General.title = 'Sin + poly(2) + Gaussian'
s.save('wobbly_peak', overwrite=True)