Model fitting¶

Specifying custom components¶

The easiest way to turn a mathematical expression into a component is using the Expression component. For example, the following is all you need to create a`Gaussian` component with more sensible parameters for spectroscopy than the one that ships with HyperSpy:

>>> g = hs.model.components1D.Expression(
... expression="height * exp(-(x - x0) ** 2 * 4 * log(2)/ fwhm ** 2)",
... name="Gaussian",
... position="x0",
... height=1,
... fwhm=1,
... centre=0,
... module="numpy")

Expression uses Sympy internally to turn the string into a funtion. By default it “translates” the expression using numpy, but often it is possible to boost performance by using numexpr instead.

Expression is only useful for analytical functions. If you know how to write the function with Python, turning it into a component is very easy modifying the following template:

from hyperspy.component import Component

class My_Component(Component):

    """
    """

    def __init__(self, parameter_1=1, parameter_2=2):
        # Define the parameters
        Component.__init__(self, ('parameter_1', 'parameter_2'))

        # Optionally we can set the initial values
         self.parameter_1.value = parameter_1
         self.parameter_1.value = parameter_1

        # The units (optional)
         self.parameter_1.units = 'Tesla'
         self.parameter_2.units = 'Kociak'

        # Once defined we can give default values to the attribute is we want
        # For example we fix the attribure_1 (optional)
         self.parameter_1.attribute_1.free = False

        # And we set the boundaries (optional)
         self.parameter_1.bmin = 0.
         self.parameter_1.bmax = None

        # Optionally, to boost the optimization speed we can define also define
        # the gradients of the function we the syntax:
        # self.parameter.grad = function
         self.parameter_1.grad = self.grad_parameter_1
         self.parameter_2.grad = self.grad_parameter_2

    # Define the function as a function of the already defined parameters, x
    # being the independent variable value
    def function(self, x):
        p1 = self.parameter_1.value
        p2 = self.parameter_2.value
        return p1 + x * p2

    # Optionally define the gradients of each parameter
     def grad_parameter_1(self, x):
         """
         Returns d(function)/d(parameter_1)
         """
         return 0

     def grad_parameter_2(self, x):
         """
         Returns d(function)/d(parameter_2)
         """
         return x

If you need help with the task please submit your question to the users mailing list.

To print the current components in a model use components of the variable. A table with component number, attribute name, component name and component type will be printed:

>>> m
<Model, title: my signal title>
>>> m.components # an empty model
   # |            Attribute Name |            Component Name |            Component Type
---- | ------------------------- | ------------------------- | -------------------------

In fact, components may be created automatically in some cases. For example, if the Signal1D is recognised as EELS data, a power-law background component will automatically be placed in the model. To add a component first we have to create an instance of the component. Once the instance has been created we can add the component to the model using the append() method, e.g. for a type of data that can be modelled using gaussians we might proceed as follows:

>>> gaussian = hs.model.components1D.Gaussian() # Create a Gaussian function component
>>> m.append(gaussian) # Add it to the model
>>> m.components # Print the model components
   # |            Attribute Name |            Component Name |            Component Type
---- | ------------------------- | ------------------------- | -------------------------
   0 |                  Gaussian |                  Gaussian |                  Gaussian
>>> gaussian2 = hs.components.Gaussian() # Create another gaussian components
>>> gaussian3 = hs.components.Gaussian() # Create a third gaussian components

We could use the append method two times to add the two gaussians, but when adding multiple components it is handier to use the extend method that enables adding a list of components at once.

>>> m.extend((gaussian2, gaussian3)) # note the double brackets!
>>> m.components
   # |            Attribute Name |            Component Name |            Component Type
---- | ------------------------- | ------------------------- | -------------------------
   0 |                  Gaussian |                  Gaussian |                  Gaussian
   1 |                Gaussian_0 |                Gaussian_0 |                  Gaussian
   2 |                Gaussian_1 |                Gaussian_1 |                  Gaussian

We can customise the name of the components.

>>> gaussian.name = 'Carbon'
>>> gaussian2.name = 'Long Hydrogen name'
>>> gaussian3.name = 'Nitrogen'
>>> m.components
   # |            Attribute Name |            Component Name |            Component Type
---- | ------------------------- | ------------------------- | -------------------------
   0 |                    Carbon |                    Carbon |                  Gaussian
   1 |        Long_Hydrogen_name |        Long Hydrogen name |                  Gaussian
   2 |                  Nitrogen |                  Nitrogen |                  Gaussian

Two components cannot have the same name.

>>> gaussian2.name = 'Carbon'
Traceback (most recent call last):
  File "<ipython-input-5-2b5669fae54a>", line 1, in <module>
    g2.name = "Carbon"
  File "/home/fjd29/Python/hyperspy/hyperspy/component.py", line 466, in name
    "the name " + str(value))
ValueError: Another component already has the name Carbon

It is possible to access the components in the model by their name or by the index in the model.

>>> m
   # |            Attribute Name |            Component Name |            Component Type
---- | ------------------------- | ------------------------- | -------------------------
   0 |                    Carbon |                    Carbon |                  Gaussian
   1 |        Long_Hydrogen_name |        Long Hydrogen name |                  Gaussian
   2 |                  Nitrogen |                  Nitrogen |                  Gaussian
>>> m[0]
<Carbon (Gaussian component)>
>>> m["Long Hydrogen name"]
<Long Hydrogen name (Gaussian component)>

In addition, the components can be accessed in the components Model attribute. This is specially useful when working in interactive data analysis with IPython because it enables tab completion.

>>> m
   # |            Attribute Name |            Component Name |            Component Type
---- | ------------------------- | ------------------------- | -------------------------
   0 |                    Carbon |                    Carbon |                  Gaussian
   1 |        Long_Hydrogen_name |        Long Hydrogen name |                  Gaussian
   2 |                  Nitrogen |                  Nitrogen |                  Gaussian
>>> m.components.Long_Hydrogen_name
<Long Hydrogen name (Gaussian component)>

It is possible to “switch off” a component by setting its active to False. When a components is switched off, to all effects it is as if it was not part of the model. To switch it on simply set the active attribute back to True.

In multidimensional signals it is possible to store the value of the active attribute at each navigation index. To enable this feature for a given component set the active_is_multidimensional attribute to True.

>>> s = hs.signals.Signal1D(np.arange(100).reshape(10,10))
>>> m = s.create_model()
>>> g1 = hs.model.components1D.Gaussian()
>>> g2 = hs.model.components1D.Gaussian()
>>> m.extend([g1,g2])
>>> g1.active_is_multidimensional = True
>>> g1._active_array
array([ True,  True,  True,  True,  True,  True,  True,  True,  True,  True], dtype=bool)
>>> g2._active_array is None
True
>>> m.set_component_active_value(False)
>>> g1._active_array
array([False, False, False, False, False, False, False, False, False, False], dtype=bool)
>>> m.set_component_active_value(True, only_current=True)
>>> g1._active_array
array([ True, False, False, False, False, False, False, False, False, False], dtype=bool)
>>> g1.active_is_multidimensional = False
>>> g1._active_array is None
True

Visualizing the model¶

To visualise the result use the plot() method:

>>> m.plot() # Visualise the results

By default only the full model line is displayed in the plot. In addition, it is possible to display the individual components by calling enable_plot_components() or directly using plot():

>>> m.plot(plot_components=True) # Visualise the results

To disable this feature call disable_plot_components().

By default the model plot is automatically updated when any parameter value changes. It is possible to suspend this feature with suspend_update(). To resume it use resume_update().

Setting the initial parameters¶

Non-linear regression often requires setting sensible starting parameters. This can be done by plotting the model and adjusting the parameters by hand.

If running in a Jupyter Notebook, interactive widgets can be used to conveniently adjust the parameter values by running notebook_interaction() for Model, Component and Parameter.

Interactive widgets for the full model in a Jupyter notebook. Drag the sliders to adjust current parameter values. Typing different minimum and maximum values changes the boundaries of the slider.

Also, enable_adjust_position() provides an interactive way of setting the position of the components with a well define position. disable_adjust_position() disables the tool.

Interactive component position adjustment tool.Drag the vertical lines to set the initial value of the position parameter.

Exclude data from the fitting process¶

The following BaseModel methods can be used to exclude undesired spectral channels from the fitting process:

• set_signal_range()
• remove_signal_range()
• reset_signal_range()

Fitting multidimensional datasets¶

To fit the model to all the elements of a multidimensional datataset use multifit(), e.g.:

>>> m.multifit() # warning: this can be a lengthy process on large datasets

multifit() fits the model at the first position, store the result of the fit internally and move to the next position until reaching the end of the dataset.

Sometimes one may like to store and fetch the value of the parameters at a given position manually. This is possible using store_current_values() and fetch_stored_values().

Visualising the result of the fit¶

The BaseModel plot_results(), Component plot() and Parameter plot() methods can be used to visualise the result of the fit when fitting multidimensional datasets.

Saving and loading the result of the fit¶

To save a model, a convenience function save() is provided, which stores the current model into its signal and saves the signal. As described in storing_models, more than just one model can be saved with one signal.

>>> m = s.create_model()
>>> # analysis and fitting goes here
>>> m.save('my_filename', 'model_name')
>>> l = hs.load('my_filename.hdf5')
>>> m = l.models.restore('model_name') # or l.models.model_name.restore()

For older versions of HyperSpy (before 0.9), the instructions were as follows:

Note that this method is known to be brittle i.e. there is no guarantee that a version of HyperSpy different from the one used to save the model will be able to load it successfully. Also, it is advisable not to use this method in combination with functions that alter the value of the parameters interactively (e.g. enable_adjust_position) as the modifications made by this functions are normally not stored in the IPython notebook or Python script.

To save a model:

1. Save the parameter arrays to a file using save_parameters2file().
2. Save all the commands that used to create the model to a file. This can be done in the form of an IPython notebook or a Python script.
3. (Optional) Comment out or delete the fitting commangs (e.g. multifit).

To recreate the model:

1. Execute the IPython notebook or Python script.
2. Use load_parameters_from_file() to load back the parameter values and arrays.

Exporting the result of the fit¶

The BaseModel export_results(), Component export() and Parameter export() methods can be used to export the result of the optimization in all supported formats.

The idea¶

The main idea of SAMFire is to change two things compared to the traditional way of fitting datasets with many dimensions in the navigation space:

1. Pick a more sensible pixel fitting order.
2. Calculate the pixel starting parameters from already fitted parts of the dataset.

Both of these aspects are linked one to another and are represented by two different strategy families that SAMFfire uses while operating.

Strategies¶

During operation SAMFire uses a list of strategies to determine how to select the next pixel and estimate its starting parameters. Only one strategy is used at a time. Next strategy is chosen when no new pixels are can be fitted with the current strategy. Once either the strategy list is exhausted or the full dataset fitted, the algorithm terminates.

There are two families of strategies. In each family there may be many strategies, using different statistical or significance measures.

As a rule of thumb, the first strategy in the list should always be from the local family, followed by a strategy from the global family.

Local strategy family¶

These strategies assume that locally neighbouring pixels are similar. As a result, the pixel fitting order seems to follow data-suggested order, and the starting values are computed from the surrounding already fitted pixels.

More information about the exact procedure will be available once the accompanying paper is published.

Global strategy family¶

Global strategies assume that the navigation coordinates of each pixel bear no relation to it’s signal (i.e. the location of pixels is meaningless). As a result, the pixels are selected at random to ensure uniform sampling of the navigation space.

A number of candidate starting values are computed form global statistical measures. These values are all attempted in order until a satisfactory result is found (not necessarily testing all available starting guesses). As a result, on average each pixel requires significantly more computations when compared to a local strategy.

More information about the exact procedure will be available once the accompanying paper is published.

Seed points¶

Due to the strategies using already fitted pixels to estimate the starting values, at least one pixel has to be fitted beforehand by the user.

The seed pixel(s) should be selected to require the most complex model present in the dataset, however in-built goodness of fit checks ensure that only sufficiently well fitted values are allowed to propagate.

If the dataset consists of regions (in the navigation space) of highly dissimilar pixels, often called “domain structures”, at least one seed pixel should be given for each unique region.

If the starting pixels were not optimal, only part of the dataset will be fitted. In such cases it is best to allow the algorithm terminate, then provide new (better) seed pixels by hand, and restart SAMFire. It will use the new seed together with the already computed parts of the data.

Usage¶

After creating a model and fitting suitable seed pixels, to fit the rest of the multi-dimensional dataset using SAMFire we must craete a SAMFire instance as follows:

>>> samf = m.create_samfire(workers=None, ipyparallel=False)

By default SAMFire will look for an ipyparallel cluster for the workers for around 30 seconds. If none is available, it will use multiprocessing instead. However, if you are not planning to use ipyparallel, it’s recommended specify it explicitly via the ipyparallel=False argument, to use the fall-back option of multiprocessing.

By default a new SAMFire object already has two (and currently only) strategies added to its strategist list:

>>> samf.strategies
  A |    # | Strategy
 -- | ---- | -------------------------
  x |    0 | Reduced chi squared strategy
    |    1 | Histogram global strategy

The currently active strategy is marked by an ‘x’ in the first column.

If a new datapoint (i.e. pixel) is added manually, the “database” of the currently active strategy has to be refreshed using the refresh_database() call.

The current strategy “database” can be plotted using the plot() method.

Whilst SAMFire is running, each pixel is checked by a goodness_test, which is by default red_chisq_test, checking the reduced chi-squared to be in the bounds of [0, 2].

This tolerance can (and most likely should!) be changed appropriately for the data as follows:

>>> samf.metadata.goodness_test.tolerance = 0.3 # use a sensible value

The SAMFire managed multi-dimensional fit can be started using the start() method. All keyword arguments are passed to the underlying (i.e. usual) fit() call:

>>> samf.start(fitter='mpfit', bounded=True)