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`

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

- General:
- For fitting
- Data -> Model
- Components
- Parameters and values
- Fitting
- After fit

- (Slightly more) Advanced:
- Generating the synthetic data:

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!

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

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

In [3]:

```
help(hs.model.components1D)
```

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]:

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]:

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]:

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]:

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

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]:

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

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

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

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

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

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

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

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

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