# 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

## Contents:¶

1. General:
2. For fitting
3. (Slightly more) Advanced:
4. Generating the synthetic data:

### Terminology and relationships¶

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¶

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?¶

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

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


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.

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


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

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

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

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


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


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¶

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


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)


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¶

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¶

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¶

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]]
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'] = 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.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¶

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

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)