MCALF: Multi-Component Atmospheric Line Fitting

MCALF Documentation

Welcome to MCALF’s documentation!

MCALF is an open-source Python package for accurately constraining velocity information from spectral imaging observations using machine learning techniques.

These pages document how the package can be interacted with. Some examples are also provided. A Documentation Index and a Module Index are available.

User Documentation

License

MCALF is an open-source Python package for accurately constraining velocity information from spectral imaging observations using machine learning techniques.

This software package is intended to be used by solar physicists trying to extract line-of-sight (LOS) Doppler velocity information from spectral imaging observations (Stokes I measurements) of the Sun. A ‘toolkit’ is provided that can be used to define a spectral model optimised for a particular dataset.

This package is particularly suited for extracting velocity information from spectral imaging observations where the individual spectra can contain multiple spectral components. Such multiple components are typically present when active solar phenomenon occur within an isolated region of the solar disk. Spectra within such a region will often have a large emission component superimposed on top of the underlying absorption spectral profile from the quiescent solar atmosphere.

A sample model is provided for an IBIS Ca II 8542 Å spectral imaging sunspot dataset. This dataset typically contains spectra with multiple atmospheric components and this package supports the isolation of the individual components such that velocity information can be constrained for each component. Using this sample model, as well as the separate base (template) model it is built upon, a custom model can easily be built for a specific dataset.

The custom model can be designed to take into account the spectral shape of each particular spectrum in the dataset. By training a neural network classifier using a sample of spectra from the dataset labelled with their spectral shapes, the spectral shape of any spectrum in the dataset can be found. The fitting algorithm can then be adjusted for each spectrum based on the particular spectral shape the neural network assigned it.

This package is designed to run in parallel over large data cubes, as well as in serial. As each spectrum is processed in isolation, this package scales very well across many processor cores. Numerous functions are provided to plot the results in a clearly. The MCALF API also contains many useful functions which have the potential of being integrated into other Python packages.

Installation

For easier package management we recommend using Miniconda (or Anaconda) and creating a new conda environment to install MCALF inside. To install MCALF using Miniconda, run the following commands in your system’s command prompt, or if you are using Windows, in the ‘Anaconda Prompt’:

$ conda config --add channels conda-forge
$ conda config --set channel_priority strict
$ conda install mcalf

MCALF is updated to the latest version by running:

$ conda update mcalf

Alternatively, you can install MCALF using pip:

$ pip install mcalf

Testing

A test suite is included with the package. The package is tested on multiple platforms, however you may wish to run the tests on your system also.

Installing Dependencies

Using MCALF with pip

To run the tests you need a number of extra packages installed. If you installed MCALF using pip, you can run pip install "mcalf[tests]" to install the additional testing dependencies (and MCALF if it’s not already installed).

Using MCALF with conda

If you want to use MCALF inside a conda environment you should first follow the conda installation instructions above. Once MCALF is installed in a conda environment, ask conda to install each of MCALF’s testing dependencies using the following command. (See setup.cfg for an up-to-date list of dependencies.)

$ conda install pytest pytest-cov pytest-mpl tox

Running Tests

Tests should be run within the virtual environment where MCALF and its testing dependencies were installed. Run the following command to test your installation,

$ pytest --pyargs mcalf

Editing the Code

If you are planning on making changes to your local version of the code, it is recommended to run the test suite to help ensure that the changes do not introduce problems elsewhere.

Before making changes, you’ll need to set up a development environment. The SunPy Community have compiled an excellent set of instructions and is available in their documentation. You can mostly replace sunpy with mcalf, and install with

$ pip install -e ".[tests,docs]"

After making changes to the MCALF source, run the MCALF test suite with the following command (while in the same directory as setup.py),

$ pytest --pyargs mcalf --cov

The tox package has also been configured to run the MCALF test suite.

Getting Started

The following examples provide the key details on how to use this package. For more details on how to use the particular classes and function, please consult the Code Reference.

Example Gallery

Here are a collection of examples on how this package can be used.

Models

Below are examples of how to use the models included within the models module:

Working with IBIS data

Working with IBIS data

Using IBIS8542Model

Using IBIS8542Model

Visualisation

Below are examples of plots produced using functions within the visualisation module:

Plot a bar chart of classifications

Plot a bar chart of classifications

Combine multiple classification plots

Combine multiple classification plots

Plot a map of classifications

Plot a map of classifications

Plot a grid of spectra grouped by classification

Plot a grid of spectra grouped by classification

Plot a fitted spectrum

Plot a fitted spectrum

Plot a map of velocities

Plot a map of velocities

Plot a spectrum

Plot a spectrum

Models

Below are examples of how to use the models included within the models module:

Working with IBIS data

Working with IBIS data

Using IBIS8542Model

Using IBIS8542Model

Working with IBIS data

This example shows how to initialise the mcalf.models.IBIS8542Model class with real IBIS data, and train a neural network classifier. We then proceed to fit the array of spectra and visualise the results.

Download sample data

First, the sample data needs to be downloaded from the GitHub repository where it is hosted. This will create four new files in the current directory (about 651 KB total). You may need to install the requests Python package for this step to run.

import requests

path = 'https://raw.githubusercontent.com/ConorMacBride/mcalf/main/examples/data/ibis8542data/'

for file in ('wavelengths.txt', 'spectra.fits',
             'training_data.json', 'results.fits'):
    r = requests.get(path + file, allow_redirects=True)
    with open(file, 'wb') as f:
        f.write(r.content)

Load the sample data

Next, the downloaded data needs to be loaded into Python.

# Import the packages needed for loading data
import json
import numpy as np
from astropy.io import fits

# Load the spectra's wavelength points
wavelengths = np.loadtxt('wavelengths.txt', dtype='>f4')

# Load the array of spectra
with fits.open('spectra.fits') as hdul:
    spectra = hdul[0].data

# Load indices of labelled spectra
with open('training_data.json', 'r') as f:
    data = f.read()
training_data = json.loads(data)

As you can see, the sample data consists of a 60 by 50 array of spectra with 27 wavelength points,

print(wavelengths.shape, spectra.shape)

(27,) (27, 60, 50)

The blue wing and line core intensity values of the spectra are plotted below for illustrative purposes,

import matplotlib.pyplot as plt
import astropy.units as u
from mcalf.visualisation import plot_map

fig, ax = plt.subplots(1, 2, sharey=True, constrained_layout=True)

wing_data = np.log(spectra[0])
core_data = np.log(spectra[len(wavelengths)//2])

res = {
    'offset': (-25, -30),
    'resolution': (0.098 * 5 * u.arcsec, 0.098 * 5 * u.arcsec),
    'show_colorbar': False,
}

wing = plot_map(wing_data, ax=ax[0], **res,
                vmin=np.min(wing_data), vmax=np.max(wing_data))
core = plot_map(core_data, ax=ax[1], **res,
                vmin=np.min(core_data), vmax=np.max(core_data))

wing.set_cmap('gray')
core.set_cmap('gray')

ax[0].set_title('blue wing')
ax[1].set_title('line core')
ax[1].set_ylabel('')

plt.show()

Generate the backgrounds array

As discussed in the Using IBIS8542Model example, a background intensity value must be specified for each spectrum.

For this small sample dataset, we shall simply use the average of the four leftmost intensity values of each spectrum,

backgrounds = np.mean(spectra[:4], axis=0)

Initialise the IBIS8542Model

The loaded data can now be passed into an mcalf.models.IBIS8542Model object.

import mcalf.models

model = mcalf.models.IBIS8542Model(original_wavelengths=wavelengths, random_state=0)

model.load_background(backgrounds, ['row', 'column'])
model.load_array(spectra, ['wavelength', 'row', 'column'])

Training the neural network

By default, the mcalf.models.IBIS8542Model object is loaded with an untrained neural network,

model.neural_network

GridSearchCV(cv=5,
             estimator=MLPClassifier(hidden_layer_sizes=(40,), max_iter=1000,
                                     random_state=0, solver='lbfgs'),
             n_jobs=-1,
             param_grid={'alpha': [1e-05, 2e-05, 3e-05, 4e-05, 5e-05, 6e-05,
                                   7e-05, 8e-05, 9e-05]})

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

The mcalf.models.IBIS8542Model class provides two methods to train and test the loaded neural network.

The training_data.json file contains a dictionary of indices for each classification 0 to 5. These indices correspond to randomly pre-chosen spectra in the spectra.fits file.

The training set consists of 200 spectra; 40 for each classification. This training set is for demonstration purposes only, generally it is not recommended to train with such a relatively high percentage of your data, as the risk of overfitting the neural network to this specific 60 by 50 array is increased.

To begin, we’ll convert the list of indices into a list of spectra and corresponding classifications,

from mcalf.utils.spec import normalise_spectrum


def select_training_set(indices, model):
    for c in sorted([int(i) for i in indices.keys()]):
        i = indices[str(c)]
        spectra = np.array([normalise_spectrum(
            model.get_spectra(row=j, column=k)[0, 0, 0],
            model.constant_wavelengths, model.constant_wavelengths
        ) for j, k in i])
        try:
            _X = np.vstack((_X, spectra))
            _y = np.hstack((_y, [c] * len(spectra)))
        except NameError:
            _X = spectra
            _y = [c] * len(spectra)
    return _X, _y


X, y = select_training_set(training_data, model)

print(X.shape)  # spectra
print(y.shape)  # labels/classifications

(200, 49)
(200,)

These classifications look as follows,

from mcalf.visualisation import plot_classifications

plot_classifications(X, y)

GridSpec(2, 3)

Now we can train the neural network on 100 labelled spectra (even indices),

model.train(X[::2], y[::2])

And now we can use the other 100 labelled spectra (odd indices) to test the performance of the neural network,

model.test(X[1::2], y[1::2])

+---------------------------------------------+
|      Neural Network Testing Statistics      |
+---------------------------------------------+
| Percentage predictions==labels ::  89.00%   |
+---------------------------------------------+
| Average deviation for each classification   |
+---------------------------------------------+
|     class  0 ::   0.00 ±  0.00              |
+---------------------------------------------+
|     class  1 ::   0.20 ±  0.68              |
+---------------------------------------------+
|     class  2 ::   0.10 ±  0.44              |
+---------------------------------------------+
|     class  3 ::   0.00 ±  0.32              |
+---------------------------------------------+
|     class  4 ::  -0.05 ±  0.22              |
+---------------------------------------------+
| Average deviation overall ::   0.05 ±  0.41 |
+---------------------------------------------+
              precision    recall  f1-score   support

           0       0.95      1.00      0.98        20
           1       0.94      0.80      0.86        20
           2       0.89      0.80      0.84        20
           3       0.75      0.90      0.82        20
           4       0.95      0.95      0.95        20

    accuracy                           0.89       100
   macro avg       0.90      0.89      0.89       100
weighted avg       0.90      0.89      0.89       100

Now that we have a trained neural network, we can use it to classify spectra. Usually spectra will be classified automatically during the fitting process, however, you can request the classification by themselves,

classifications = model.classify_spectra(row=range(60), column=range(50))

These classifications look as follows,

from mcalf.visualisation import plot_class_map

plot_class_map(classifications)

<matplotlib.image.AxesImage object at 0x7f0cc7869580>

Creating a reproducible classifier

The neural network classifier introduces a certain amount of randomness when it it fitting based on the training data. This randomness arises in the initial values of the weights and biases that are fitted during the training process, as well as the order in which the training data are used.

This means that two neural networks trained on identical data will not produce the same results. To aid the reproducibility of results that rely on a neural network’s classifications, a random_state integer can be passed to mcalf.models.IBIS8542Model as we did above. When we set this value to an integer, no matter how many times we train the neural network on the same data, it will always give the same results.

Until better solutions are available to store trained neural networks, a trained neural network can be saved to a Python pickle file and later reloaded. For maximum compatibility, it is recommended to reload into the same version of scikit-learn and its dependencies.

The neural network trained above can be saved as follows,

import pickle
pkl = open('trained_neural_network.pkl', 'wb')
pickle.dump(model.neural_network, pkl)
pkl.close()

This trained neural network can then be reloaded at a later date as follows,

import pickle
pkl = open('trained_neural_network.pkl', 'rb')
model.neural_network = pickle.load(pkl)  # Overwrite the default untrained model

And you can see that the classifications of spectra are the same,

plot_class_map(model.classify_spectra(row=range(60), column=range(50)))

<matplotlib.image.AxesImage object at 0x7f0cc6db5280>

Please see the scikit-learn documentation for more details on model persistence.

Fitting the spectra

Now that the data have been loaded and the neural network has been trained, we can proceed to fit the spectra.

Using pre-calculated results

As our 60 by 50 array contains 3000 spectra, it would take roughly 10 minutes to fit them all over 6 processing pools. We include pre-calculated results in the downloaded results.fits file.

The next step of the example loads this file back into Python as though we have just directly calculated it. This isn’t something you would usually need to do, so do not worry about the contents of the load_results() function, however, we plan to include this functionality in MCALF itself in the future.

def load_results(file):

    with fits.open(file) as hdul:
        for hdu in hdul:
            if hdu.name == 'PARAMETERS':
                r_parameters = hdu.data.copy().reshape(-1, 8)
            elif hdu.name == 'CLASSIFICATIONS':
                r_classifications = hdu.data.copy().flatten()
                r_profile = np.full_like(r_classifications, 'both', dtype=object)
                r_profile[r_classifications <= 1] = 'absorption'
            elif hdu.name == 'SUCCESS':
                r_success = hdu.data.copy().flatten()
            elif hdu.name == 'CHI2':
                r_chi2 = hdu.data.copy().flatten()

    results = []
    for i in range(len(r_parameters)):
        fitted_parameters = r_parameters[i]
        fit_info = {
            'classification': r_classifications[i],
            'profile': r_profile[i],
            'success': r_success[i],
            'chi2': r_chi2[i],
            'index': [0, *np.unravel_index(i, (60, 50))],
        }
        if fit_info['profile'] == 'absorption':
            fitted_parameters = fitted_parameters[:4]
        results.append(mcalf.models.FitResult(fitted_parameters, fit_info))

    return results


result_list = load_results('results.fits')

result_list[:4]  # The first four

[Successful FitResult at (0, 0, 0) with absorption profile of classification 0, Successful FitResult at (0, 0, 1) with absorption profile of classification 0, Successful FitResult at (0, 0, 2) with absorption profile of classification 0, Successful FitResult at (0, 0, 3) with absorption profile of classification 0]

Using your own results

You can run the following code to generate the result_list variable for yourself. Try starting with a smaller range of rows and columns, and set the number of pools based on the specification of your processor.

# result_list = model.fit(row=range(60), column=range(50), n_pools=6)

The order of the mcalf.models.FitResult objects in this list will also differ as the order that spectra finish fitting in each pool is unpredictable.

Merging the FitResult objects

The list of mcalf.models.FitResult objects can be mergerd into a mcalf.models.FitResults object, and then saved to a file, just like the results.fits file downloaded earlier.

First the object needs to be initialised with the spatial dimensions and number of fitted parameters,

results = mcalf.models.FitResults((60, 50), 8)

Now we can loop through the list of fits and append them to this object,

for fit in result_list:
    results.append(fit)

This object can then be saved to file,

results.save('ibis8542data.fits', model)

The file has the following structure,

with fits.open('ibis8542data.fits') as hdul:
    hdul.info()

Filename: ibis8542data.fits
No.    Name      Ver    Type      Cards   Dimensions   Format
  0  PRIMARY       1 PrimaryHDU      13   (0,)
  1  PARAMETERS    1 ImageHDU        12   (8, 50, 60)   float64
  2  CLASSIFICATIONS    1 ImageHDU        10   (50, 60)   int16
  3  PROFILE       1 ImageHDU        11   (50, 60)   int16
  4  SUCCESS       1 ImageHDU        10   (50, 60)   int16
  5  CHI2          1 ImageHDU        10   (50, 60)   float64
  6  VLOSA         1 ImageHDU        12   (50, 60)   float64
  7  VLOSQ         1 ImageHDU        12   (50, 60)   float64

Exploring the fitted data

You can plot a fitted spectrum as follows,

model.plot(result_list[0])

<Axes: xlabel='wavelength (Å)', ylabel='intensity'>

You can calculate and plot Doppler velocities for both the quiescent and active regimes as follows, (with an outline of the sunspot’s umbra),

vq = results.velocities(model, vtype='quiescent')
va = results.velocities(model, vtype='active')

umbra_mask = np.full_like(backgrounds, True)
umbra_mask[backgrounds < 1100] = False

fig, ax = plt.subplots(1, 2, sharey=True, constrained_layout=True)

settings = {
    'show_colorbar': False, 'vmax': 4, 'offset': (-25, -30),
    'resolution': (0.098 * 5 * u.arcsec, 0.098 * 5 * u.arcsec),
    'umbra_mask': umbra_mask,
}

im = plot_map(vq, ax=ax[0], **settings)
plot_map(va, ax=ax[1], **settings)

ax[0].set_title('quiescent')
ax[1].set_title('active')
ax[1].set_ylabel('')
fig.colorbar(im, ax=ax, location='bottom', label='Doppler velocity (km/s)')

plt.show()

Using IBIS8542Model

This is an example showing how to fit an array of spectra using the mcalf.models.IBIS8542Model class, and plot and export the results.

Generate sample data

Create spectral grid

First, we shall generate some random sample data to demonstrate the API. The randomly generated spectra are not intended to be representative of observed spectra numerically, they just have a similar shape.

import numpy as np
import matplotlib.pyplot as plt

from mcalf.profiles.voigt import voigt, double_voigt
from mcalf.visualisation import plot_spectrum

# Create demo wavelength grid
w = np.linspace(8541, 8543, 20)

We shall generate demo spectra in a 2x3 grid. Half of the spectra will have one spectral component (absorption only) and the other half will have two spectral components (mixed absorption and emission components).

Demo spectral components will be modelled as Voigt functions with randomly generated parameters, each within a set range of values.

def v(classification, w, *args):
    """Voigt function wrapper."""
    if classification < 1.5:  # absorption only
        return voigt(w, *args[:4], args[-1])
    return double_voigt(w, *args)  # absorption + emission


def s():
    """Generate random spectra for a 2x3 grid."""
    np.random.seed(0)  # same spectra every time
    p = np.random.rand(9, 6)  # random Voigt parameters
    # 0  1  2  3  4  5  6  7  8  # p index
    # a1 b1 s1 g1 a2 b2 s2 g2 d  # Voigt parameter
    # absorption |emission   |background

    p[0] = 100 * p[0] - 1000  # absorption amplitude
    p[4] = 100 * p[4] + 1000  # emission amplitude

    for i in (1, 5):  # abs. and emi. peak positions
        p[i] = 0.05 * p[i] - 0.025 + 8542

    for i in (2, 3, 6, 7):  # Voigt sigmas and gammas
        p[i] = 0.1 * p[i] + 0.1

    p[8] = 300 * p[8] + 2000  # intensity background constant

    # Define each spectrum's classification
    c = [0, 2, 0, 2, 0, 2]
    # Choose single or double component spectrum
    # based on this inside the function `v()`.

    # Generate the spectra
    specs = [v(c[i], w, *p[:, i]) for i in range(6)]

    # Reshape to 2x3 grid
    return np.asarray(specs).reshape((2, 3, len(w)))


raw_data = s()

print('shape of spectral grid:', raw_data.shape)

shape of spectral grid: (2, 3, 20)

These spectra look as follows,

fig, axes = plt.subplots(2, 3, constrained_layout=True)
for ax, spec in zip(axes.flat, raw_data.reshape((6, raw_data.shape[-1]))):

    plot_spectrum(w, spec, normalised=False, ax=ax)

plt.show()

Compute the background intensity

MCALF does not model a constant background value, i.e., the fitting process assumes the intensity values far out in the spectrum’s wings are zero.

You can however tell MCALF what the constant background value is, and it will subtract it from the spectrum before fitting.

This process was not made automatic as we wanted to give the user full control on setting the background value.

For these demo data, we shall simply set the background to the first intensity value of every spectrum. For a real dataset, you may wish to average the background value throughout a range of spectral points, or even do a moving average throughout time. Functions are provided in mcalf.utils.smooth to assist with this process.

backgrounds = raw_data[:, :, 0]

print('shape of background intensity grid:', backgrounds.shape)

shape of background intensity grid: (2, 3)

Define a demo classifier

In this example we will not demonstrate how to create a neural network classifier. MCALF offers a lot of flexibility when it comes to the classifier. This demo classifier outlines the basic API that is required. By default, the model has a sklearn.neural_network.MLPClassifier object preloaded for use as the classifier. The methods mcalf.models.ModelBase.train() and mcalf.models.ModelBase.test() are provided by MCALF to assist with training the neural network. There is also a useful script in the Getting Started section under User Documentation in the sidebar for semi-automating the process of creating the ground truth dataset.

Please see tutorials for packages such as scikit-learn for more in-depth advice on creating classifiers.

As we only have six spectra with two distinct shapes, we can create a very simple classifier that classifies spectra based on whether their central intensity is greater or smaller than the left most intensity.

class DemoClassifier:
    @staticmethod
    def predict(X):
        y = np.zeros(len(X), dtype=int)
        y[X[:, len(X[0]) // 2] > X[:, 0]] = 2
        return y

Using IBIS8542Model with the generated data

Initialise the model

Everything we have been doing up to this point has been creating the demo data and classifier. Now we can actually create a model and load in the demo data, although the following steps would be identical for a real dataset.

from mcalf.models import IBIS8542Model

# Initialise the model with the wavelength grid
model = IBIS8542Model(original_wavelengths=w)

# Load the spectral shape classifier
model.neural_network = DemoClassifier

# Load the array of spectra and background intensities
model.load_array(raw_data, ['row', 'column', 'wavelength'])
model.load_background(backgrounds, ['row', 'column'])

Fit the loaded data

We have now fully initialised the model. We can now call methods to fit the model to the loaded spectra. In the following step we fit all the loaded spectra and a 1D list of FitResult objects is returned.

fits = model.fit(row=[0, 1], column=range(3))

print(fits)

Processing 6 spectra
/home/docs/checkouts/readthedocs.org/user_builds/mcalf/checkouts/stable/src/mcalf/models/base.py:1003: UserWarning: RuntimeError(Optimal parameters not found: The maximum number of function evaluations is exceeded.) at (0, 0, 1)
  warnings.warn(f"RuntimeError({e}){location_text}")
/home/docs/checkouts/readthedocs.org/user_builds/mcalf/checkouts/stable/src/mcalf/models/base.py:1003: UserWarning: RuntimeError(Optimal parameters not found: The maximum number of function evaluations is exceeded.) at (0, 1, 2)
  warnings.warn(f"RuntimeError({e}){location_text}")
[Successful FitResult with absorption profile of classification 0, Unsuccessful FitResult with both profile of classification 2, Successful FitResult with absorption profile of classification 0, Successful FitResult with both profile of classification 2, Successful FitResult with absorption profile of classification 0, Unsuccessful FitResult with both profile of classification 2]

Plot the fitted parameters

The individual components of each fit can now be plotted separately,

fig, axes = plt.subplots(2, 3, constrained_layout=True)
for ax, fit in zip(axes.flat, fits):

    model.plot_separate(fit, show_legend=False, ax=ax)

plt.show()

Extract processed spectra

As well as fitting spectra, we can call other methods. In this step we’ll extract the array of loaded spectra. However, these spectra have been re-interpolated to a new finer wavelength grid. (This grid can be customised when initialising the model.)

spectra = model.get_spectra(row=[0, 1], column=range(3))

print('new shape of spectral grid:', spectra.shape)

spectra_1d = spectra.reshape((6, spectra.shape[-1]))

fig, axes = plt.subplots(2, 3, constrained_layout=True)
for ax, spec in zip(axes.flat, spectra_1d):

    plot_spectrum(model.constant_wavelengths, spec,
                  normalised=False, ax=ax)

plt.show()

new shape of spectral grid: (1, 2, 3, 41)

Classify spectra

We can also classify the loaded spectra and create plots,

classifications = model.classify_spectra(row=[0, 1], column=range(3))

print('classifications are:', classifications)

classifications are: [[[0 2 0]
  [2 0 2]]]

This function plots the spectra grouped by classification,

from mcalf.visualisation import plot_classifications

plot_classifications(spectra_1d, classifications.flatten())

GridSpec(1, 2)

This function plots a spatial map of the classifications on the 2x3 grid,

from mcalf.visualisation import plot_class_map

plot_class_map(classifications)

<matplotlib.image.AxesImage object at 0x7f0cc53bb6d0>

Merge output data

The FitResult objects in the fits 1D list can be merged into a grid. Each of these objects can be appended to a single FitResults object. This allows for increased data portability,

from mcalf.models import FitResults

# Initialise with the grid shape and num. params.
results = FitResults((2, 3), 8)

# Append each fit to the object
for fit in fits:
    results.append(fit)

Now that we have appended all 6 fits, we can access the merged data,

print(results.classifications)
print(results.profile)
print(results.success)
print(results.parameters)

[[0 2 0]
 [2 0 2]]
[['absorption' 'both' 'absorption']
 ['both' 'absorption' 'both']]
[[ True False  True]
 [ True  True False]]
[[[-8.12854071e+02  8.54199697e+03  1.78499712e-01  1.26560705e-01
               nan             nan             nan             nan]
  [            nan             nan             nan             nan
               nan             nan             nan             nan]
  [-8.18417945e+02  8.54202321e+03  1.28952238e-01  1.56493579e-01
               nan             nan             nan             nan]]

 [[-2.21478736e+03  8.54199483e+03  1.06846047e-01  2.21738766e-01
    2.33701044e+03  8.54199968e+03  1.48605642e-01  2.05127799e-01]
  [-8.73401645e+02  8.54201460e+03  1.15790594e-01  1.20586703e-01
               nan             nan             nan             nan]
  [            nan             nan             nan             nan
               nan             nan             nan             nan]]]

Calculate velocities

And finally, we can calculate Doppler velocities for both the quiescent (absorption) and active (emission) regimes. (The model needs to be given so the stationary line core wavelength is available.)

quiescent = results.velocities(model, vtype='quiescent')
active = results.velocities(model, vtype='active')

from mcalf.visualisation import plot_map

fig, axes = plt.subplots(2, constrained_layout=True)
plot_map(quiescent, ax=axes[0])
plot_map(active, ax=axes[1])
plt.show()

Export to FITS file

The FitResults object can also be exported to a FITS file,

results.save('ibis8542model_demo_fit.fits', model)

The file has the following structure,

from astropy.io import fits

with fits.open('ibis8542model_demo_fit.fits') as hdul:
    hdul.info()

Filename: ibis8542model_demo_fit.fits
No.    Name      Ver    Type      Cards   Dimensions   Format
  0  PRIMARY       1 PrimaryHDU      13   (0,)
  1  PARAMETERS    1 ImageHDU        12   (8, 3, 2)   float64
  2  CLASSIFICATIONS    1 ImageHDU        10   (3, 2)   int16
  3  PROFILE       1 ImageHDU        11   (3, 2)   int16
  4  SUCCESS       1 ImageHDU        10   (3, 2)   int16
  5  CHI2          1 ImageHDU        10   (3, 2)   float64
  6  VLOSA         1 ImageHDU        12   (3, 2)   float64
  7  VLOSQ         1 ImageHDU        12   (3, 2)   float64

This example outlines the basics of using the mcalf.models.IBIS8542Model class. Typically millions of spectra would need to be fitted and processed at the same time. For more details on running big jobs and how spectra can be fitted in parallel, see the Working with IBIS data example in the gallery.

Note

Due to limitations with the documentation hosting provider, the examples in the MCALF documentation are computed using a Python based implementation of the Voigt profile instead of the more efficient C version. When this code is run on your own machine it should take much less time than the value reported below.

Visualisation

Below are examples of plots produced using functions within the visualisation module:

Plot a bar chart of classifications

Plot a bar chart of classifications

Combine multiple classification plots

Combine multiple classification plots

Plot a map of classifications

Plot a map of classifications

Plot a grid of spectra grouped by classification

Plot a grid of spectra grouped by classification

Plot a fitted spectrum

Plot a fitted spectrum

Plot a map of velocities

Plot a map of velocities

Plot a spectrum

Plot a spectrum

Plot a bar chart of classifications

This is an example showing how to produce a bar chart showing the percentage abundance of each classification in a 2D or 3D array of classifications.

First we shall create a random 3D grid of classifications that can be plotted. Usually you would use a method such as mcalf.models.ModelBase.classify_spectra() to classify an array of spectra.

from mcalf.tests.helpers import class_map as c

t = 3  # Three images
x = 50  # 50 coordinates along x-axis
y = 20  # 20 coordinates along y-axis
n = 5  # Possible classifications [0, 1, 2, 3, 4]

class_map = c(t, x, y, n)  # 3D array of classifications (t, y, x)

Next, we shall import mcalf.visualisation.bar().

from mcalf.visualisation import bar

We can now simply plot the 3D array. By default, the first dimension of a 3D array will be averaged to produce a time average, selecting the most common classification at each (x, y) coordinate. This means the percentage abundances will correspond to the most common classification at each coordinate.

bar(class_map)

<BarContainer object of 5 artists>

Instead, the percentage abundances can be determined for the whole 3D array of classifications by setting reduce=True. This skips the averaging process.

bar(class_map, reduce=False)

<BarContainer object of 5 artists>

Alternatively, a 2D array can be passed to the function. If a 2D array is passed, no averaging is needed, and the reduce parameter is ignored.

A narrower range of classifications to be plotted can be requested with the vmin and vmax parameters. To show bars for only classifcations 1, 2, and 3,

bar(class_map, vmin=1, vmax=3)

<BarContainer object of 3 artists>

An alternative set of colours can be requested. Passing a name of a matplotlib colormap to the style parameter will produce a corresponding list of colours for each of the bars. For advanced use, explore the cmap parameter.

bar(class_map, style='viridis')

<BarContainer object of 5 artists>

The bar function integrates well with matplotlib, allowing extensive flexibility.

import matplotlib.pyplot as plt

fig, ax = plt.subplots(1, 2, constrained_layout=True)

bar(class_map, vmax=2, style='viridis', ax=ax[0])
bar(class_map, vmin=3, style='cividis', ax=ax[1])

ax[0].set_title('first 3')
ax[1].set_title('last 2')

ax[1].set_ylabel('')

plt.show()

Note that vmin and vmax are applied before the reduce value is applied. So setting these ranges can change the calculated abundances for other classifications if class_map is 3D and reduce=True.

The bars do not add up to 100% as a bar for negative, invalid classifications (and therefore classifications out of the vmin and vmax range) is not shown.

Combine multiple classification plots

This is an example showing how to use multiple classification plotting functions in a single figure.

First we shall create a random 3D grid of classifications that can be plotted. Usually you would use a method such as mcalf.models.ModelBase.classify_spectra() to classify an array of spectra.

from mcalf.tests.helpers import class_map as c

t = 1  # One image
x = 20  # 20 coordinates along x-axis
y = 20  # 20 coordinates along y-axis
n = 3  # Possible classifications [0, 1, 2]

class_map = c(t, x, y, n)  # 3D array of classifications (t, y, x)

Next we shall create a random array of spectra each labelled with a random classifications. Usually you would provide your own set of hand labelled spectra taken from spectral imaging observations of the Sun. Or you could provide a set of spectra labelled by the classifier.

from mcalf.tests.visualisation.test_classifications import spectra as s

n = 400  # 200 spectra
w = 20  # 20 wavelength points for each spectrum
low, high = 0, 3  # Possible classifications [0, 1, 2]

# 2D array of spectra (n, w), 1D array of labels (n,)
spectra, labels = s(n, w, low, high)

If a GridSpec returned by the plot_classification function has free space, a new axes can be added to the returned GridSpec. We can then request plot_class_map to plot onto this new axes. The colorbar axes can be set to fig.axes such that the colorbar takes the full height of the figure, as in this case, its colours are the same as the line plots.

import matplotlib.pyplot as plt
from mcalf.visualisation import plot_classifications, plot_class_map

fig = plt.figure(constrained_layout=True)

gs = plot_classifications(spectra, labels, nrows=2, show_labels=False)

ax = fig.add_subplot(gs[-1])
plot_class_map(class_map, ax=ax, colorbar_settings={
    'ax': fig.axes,
    'label': 'classification',
})

<matplotlib.image.AxesImage object at 0x7f0cc6ac1fd0>

The function mcalf.visualisation.init_class_data`() is intended to be an internal function for generating data that is common to multiple plotting functions. However, it may be used externally if necessary.

from mcalf.visualisation import init_class_data, bar

fig, ax = plt.subplots(1, 2, constrained_layout=True)

data = init_class_data(class_map, resolution=(0.25, 0.25), ax=ax[1])

bar(data=data, ax=ax[0])
plot_class_map(data=data, ax=ax[1], show_colorbar=False)

<matplotlib.image.AxesImage object at 0x7f0cc693eca0>

The following example should be equivalent to the example above,

fig, ax = plt.subplots(1, 2, constrained_layout=True)

bar(class_map, ax=ax[0])
plot_class_map(class_map, ax=ax[1], show_colorbar=False,
               resolution=(0.25, 0.25))

<matplotlib.image.AxesImage object at 0x7f0cc51b41c0>

Plot a map of classifications

This is an example showing how to produce a map showing the spatial distribution of spectral classifications in a 2D region of the Sun.

First we shall create a random 3D grid of classifications that can be plotted. Usually you would use a method such as mcalf.models.ModelBase.classify_spectra() to classify an array of spectra.

from mcalf.tests.helpers import class_map as c

t = 3  # Three images
x = 50  # 50 coordinates along x-axis
y = 20  # 20 coordinates along y-axis
n = 5  # Possible classifications [0, 1, 2, 3, 4]

class_map = c(t, x, y, n)  # 3D array of classifications (t, y, x)

Next, we shall import mcalf.visualisation.plot_class_map().

from mcalf.visualisation import plot_class_map

We can now simply plot the 3D array. By default, the first dimension of a 3D array will be averaged to produce a time average, selecting the most common classification at each (x, y) coordinate.

plot_class_map(class_map)

<matplotlib.image.AxesImage object at 0x7f0cc568d100>

A spatial resolution with units can be specified for each axis.

import astropy.units as u

plot_class_map(class_map, resolution=(0.75 * u.km, 1.75 * u.Mm),
               offset=(-25, -10),
               dimension=('x distance', 'y distance'))

<matplotlib.image.AxesImage object at 0x7f0cc7cb3640>

A narrower range of classifications to be plotted can be requested with the vmin and vmax parameters. Classifications outside of the range will appear as grey, the same as pixels with a negative, unassigned classification.

plot_class_map(class_map, vmin=1, vmax=3)

<matplotlib.image.AxesImage object at 0x7f0cc56135b0>

An alternative set of colours can be requested. Passing a name of a matplotlib colormap to the style parameter will produce a corresponding list of colours for each of the classifications. For advanced use, explore the cmap parameter.

plot_class_map(class_map, style='viridis')

<matplotlib.image.AxesImage object at 0x7f0cc6d361f0>

The plot_class_map function integrates well with matplotlib, allowing extensive flexibility. This example also shows how you can plot a 2D class_map and skip the averaging.

import matplotlib.pyplot as plt

fig, ax = plt.subplots(2, constrained_layout=True)

plot_class_map(class_map[0], style='viridis', ax=ax[0],
               show_colorbar=False)
plot_class_map(class_map[1], style='viridis', ax=ax[1],
               colorbar_settings={'ax': ax, 'label': 'classification'})

ax[0].set_title('time $t=0$')
ax[1].set_title('time $t=1$')

plt.show()

Plot a grid of spectra grouped by classification

This is an example showing how to produce a grid of line plots of an array of spectra labelled with a classification.

First we shall create a random array of spectra each labelled with a random classifications. Usually you would provide your own set of hand labelled spectra taken from spectral imaging observations of the Sun.

from mcalf.tests.visualisation.test_classifications import spectra as s

n = 200  # 200 spectra
w = 20  # 20 wavelength points for each spectrum
low, high = 1, 5  # Possible classifications [1, 2, 3, 4]

# 2D array of spectra (n, w), 1D array of labels (n,)
spectra, labels = s(n, w, low, high)

Next, we shall import mcalf.visualisation.plot_classifications().

from mcalf.visualisation import plot_classifications

We can now plot a simple grid of the spectra grouped by their classification. By default, a maximum of 20 spectra are plotted for each classification. The first 20 spectra are selected.

plot_classifications(spectra, labels)

GridSpec(2, 2)

A specific number of rows or columns can be requested,

plot_classifications(spectra, labels, ncols=1)

GridSpec(4, 1)

The plot settings can be adjusted,

plot_classifications(spectra, labels, show_labels=False, nlines=5,
                     style='viridis', plot_settings={'ls': '--', 'marker': 'x'})

GridSpec(2, 2)

By default, the y-axis goes from 0 to 1. This is because labelled training data will typically be rescaled between 0 and 1. However, if a particular classification has spectra that are not within 0 and 1, the y-axis limits are determined by matplotlib.

spectra[labels == 2] += 0.5
plot_classifications(spectra, labels)

GridSpec(2, 2)

Plot a fitted spectrum

This is an example showing how to plot the result of fitting a spectrum using the IBIS8542Model class.

First we shall create a list of wavelengths, with a variable wavelength spacing. Next, we shall use the Voigt profile to generate spectral intensities at each of the wavelength points. Typically you would provide a spectrum obtained from observations.

The data in this example are produced from randomly selected parameters, so numerical values in this example should be ignored.

import numpy as np
from mcalf.models import IBIS8542Model
from mcalf.profiles.voigt import double_voigt

# Create the wavelength grid and intensity values
wavelengths = np.linspace(8541, 8543, 20)
wavelengths = np.delete(wavelengths, np.s_[1:6:2])
wavelengths = np.delete(wavelengths, np.s_[-6::2])

spectrum = double_voigt(wavelengths, -526, 8542, 0.1, 0.1,
                        300, 8541.9, 0.2, 0.05, 1242)

from mcalf.visualisation import plot_spectrum

plot_spectrum(wavelengths, spectrum, normalised=False)

<Axes: xlabel='wavelength (Å)', ylabel='intensity ($I$)'>

A basic model is created,

model = IBIS8542Model(original_wavelengths=wavelengths)

The spectrum is provided to the model’s fit method. A classifier has not been loaded so the classification must be provided manually. The fitting algorithm assumes that the intensity at the ends of the spectrum is zero, so in this case we need to provide it with a background value to subtract from the spectrum before fitting.

fit = model.fit(spectrum=spectrum, classifications=4, background=1242)

print(fit)

Successful FitResult with both profile of classification 4

The spectrum can now be plotted,

model.plot(fit, spectrum=spectrum, background=1242)

<Axes: xlabel='wavelength (Å)', ylabel='intensity'>

If an array of spectra and associated background values had been loaded into the model with the load_array() and load_background() methods respectively, the spectrum and background parameters would not have to be specified when plotting. This is because the fit object would contain indices that the model object would use to look up the original loaded values.

Equivalent to above, the plot method can be called on the fit object directly. Remember to specify the model which is needed for additional information such as the stationary line core value.

fit.plot(model, spectrum=spectrum, background=1242)

If the fit has multiple spectral components, such as an active emission profile mixed with a quiescent absorption profile, the follow method can be used to plot the components separatly.

If the fit only has a single component the plot method as shown above is used.

model.plot_separate(fit, spectrum=spectrum, background=1242)

If the fit has an emission component, it is subtracted from the raw spectral data. Otherwise, the default plot method is used.

model.plot_subtraction(fit, spectrum=spectrum, background=1242)

The same line on multiple plots is only labelled the first time it plotted in the figure. This prevents duplicated entries in legends.

import matplotlib.pyplot as plt

fig, ax = plt.subplots(1, 2)

model.plot(fit, spectrum=spectrum, background=1242,
           show_legend=False, ax=ax[0])
model.plot_separate(fit, spectrum=spectrum, background=1242,
                    show_legend=False, ax=ax[1])

fig.subplots_adjust(top=0.75)  # Create space above for legend
fig.legend(ncol=2, loc='upper center', bbox_to_anchor=(0.5, 0.97))

plt.show()

The underlying mcalf.visualisation.plot_ibis8542() function can be used directly. However, it is recommended to plot using the method detailed above as it will do additional processing to the wavelengths and spectrum and also pass additional parameters, such as sigma, to this fitting function.

from mcalf.visualisation import plot_ibis8542

plot_ibis8542(wavelengths, spectrum, fit.parameters, 1242)

<Axes: xlabel='wavelength (Å)', ylabel='intensity'>

The y-axis and legend can be easily hidden,

model.plot(fit, spectrum=spectrum, background=1242,
           show_intensity=False, show_legend=False)

<Axes: xlabel='wavelength (Å)'>

Plot a map of velocities

This is an example showing how to produce a map showing the spatial distribution of velocities in a 2D region of the Sun.

First we shall create a random 2D grid of velocities that can be plotted. Usually you would use a method such as mcalf.models.FitResults.velocities() to extract an array of velocities from fitted spectra.

import numpy as np
np.random.seed(0)

x = 50  # 50 coordinates along x-axis
y = 40  # 40 coordinates along y-axis
low, high = -10, 10  # range of velocities (km/s)


def a(x, y, low, high):
    arr = np.random.normal(0, (high - low) / 2 * 0.3, (y, x))
    arr[arr < low] = low
    arr[arr > high] = high
    i = np.random.randint(0, arr.size, arr.size // 100)
    arr[np.unravel_index(i, arr.shape)] = np.nan
    return arr


arr = a(x, y, low, high)  # 2D array of velocities (y, x)

Next, we shall import mcalf.visualisation.plot_map().

from mcalf.visualisation import plot_map

We can now simply plot the 2D array.

plot_map(arr)

<matplotlib.image.AxesImage object at 0x7f0cc79778e0>

Notice that pixels with missing data (NaN) are shown in grey.

By default, the velocity data are assumed to have units km/s. If your data are not in km/s, you must either 1) rescale the array such that it is in km/s, 2) attach an astropy unit to the array to override the default, or 3) pass an astropy unit to the unit parameter to override the default. For example, we can change from km/s to m/s,

import astropy.units as u

plot_map(arr * 1000 * u.m / u.s)

<matplotlib.image.AxesImage object at 0x7f0cc7ac0dc0>

A spatial resolution with units can be specified for each axis.

plot_map(arr, resolution=(0.5 * u.km, 0.6 * u.Mm), offset=(-25, -20))

<matplotlib.image.AxesImage object at 0x7f0cc6cfc820>

A narrower range of velocities to be plotted can be requested with the vmin and vmax parameters. Classifications outside of the range will appear saturated. Providing only one of vmin and vmax with set the other such that zero is the midpoint.

plot_map(arr, vmax=4)

<matplotlib.image.AxesImage object at 0x7f0cc79f0e20>

A mask can be applied to the velocity array to isolate a region of interest. This functionally is useful if, for example, data only exist for a circular region and you want to distinguish between the pixels that are out of bounds and the data that were not successfully fitted.

from mcalf.utils.mask import genmask
mask = genmask(50, 40, 18)

plot_map(arr, ~mask)

<matplotlib.image.AxesImage object at 0x7f0cc69931c0>

Notice how data out of bounds are grey, while data which were not fitted successfully are now black.

A region of interest, typically the umbra of a sunspot, can be outlined by passing a different mask.

umbra_mask = genmask(50, 40, 5, 5, 5)

plot_map(arr, ~mask, umbra_mask)

<matplotlib.image.AxesImage object at 0x7f0cc5522bb0>

The plot_map function integrates well with matplotlib, allowing extensive flexibility.

import matplotlib.pyplot as plt

fig, ax = plt.subplots(1, 2, constrained_layout=True)

plot_map(arr[:, :25], vmax=6, ax=ax[0], show_colorbar=False)
im = plot_map(arr[:, 25:], vmax=6, ax=ax[1], show_colorbar=False)

fig.colorbar(im, ax=ax, location='bottom', label='velocity (km/s)')

ax[0].set_title('first half')
ax[1].set_title('second half')

plt.show()

Plot a spectrum

This is an example showing how to plot a spectrum with the mcalf.visualisation.plot_spectrum() function.

First we shall create a list of wavelengths, with a variable wavelength spacing. Next, we shall use the Voigt profile to generate spectral intensities at each of the wavelength points. Typically you would provide a spectrum obtained from observations.

import numpy as np
wavelengths = np.linspace(8541, 8543, 20)
wavelengths = np.delete(wavelengths, np.s_[1:6:2])
wavelengths = np.delete(wavelengths, np.s_[-6::2])

from mcalf.profiles.voigt import voigt
spectrum = voigt(wavelengths, -526, 8542, 0.1, 0.1, 1242)

Next, we shall import mcalf.visualisation.plot_spectrum().

from mcalf.visualisation import plot_spectrum

We can now simply plot the spectrum.

plot_spectrum(wavelengths, spectrum)

<Axes: xlabel='wavelength (Å)', ylabel='normalised intensity ($I/I_c$)'>

Notice how the spectrum above is normalised. The normalisation is applied by dividing through by the mean of the three rightmost points. To plot the raw spectrum,

plot_spectrum(wavelengths, spectrum, normalised=False)

<Axes: xlabel='wavelength (Å)', ylabel='intensity ($I$)'>

The line connecting the points provided in the spectrum array above is smooth. This is due to spline interpolation being applied. Interpolation can be disabled, resulting in a straight line between each of the points.

plot_spectrum(wavelengths, spectrum, smooth=False)

<Axes: xlabel='wavelength (Å)', ylabel='normalised intensity ($I/I_c$)'>

example1: Basic usage of the package

FittingIBIS.ipynb

• View code

• Download FittingIBIS.ipynb

This file is an IPython Notebook containing examples of how to use the package to accomplish typical tasks.

FittingIBIS.pro

• Download FittingIBIS.pro

This file is similar to FittingIBIS.ipynb file, except it written is IDL. It is not recommended to use the IDL wrapper in production, just use it to explore the code if you are familiar with IDL and not Python. If you wish to use this package, please use the Python implementation. IDL is not fully supported in the current version of the code for reasons such as, the Python tuple datatype cannot be passed from IDL to Python, resulting in certain function calls not being possible.

config.yml

• Download config.yml

This is an example configuration file containing default parameters. This can be easier than setting the parameters in the code. The file follows the YAML format.

Labelling Tutorial

This Jupyter notebook provides a simple, semi-automated, method to produce a ground truth data set that can be used to train a neural network for use as a spectral shape classifier in the MCALF package. The following code can be adapted depending on the number of classifications that you want.

Download LabellingTutorial.ipynb

Load the required packages

import mcalf.models
from mcalf.utils.spec import normalise_spectrum
import requests
import numpy as np
from astropy.io import fits
import matplotlib.pyplot as plt

Download sample data

path = 'https://raw.githubusercontent.com/ConorMacBride/mcalf/main/examples/data/ibis8542data/'

for file in ('wavelengths.txt', 'spectra.fits'):
    r = requests.get(path + file, allow_redirects=True)
    with open(file, 'wb') as f:
        f.write(r.content)

Load data files

wavelengths = np.loadtxt('wavelengths.txt')  # Original wavelengths

with fits.open('spectra.fits') as hdul:  # Raw spectral data
    datacube = np.asarray(hdul[0].data, dtype=np.float64)

Initialise the model that will use the labelled data

model = mcalf.models.IBIS8542Model(original_wavelengths=wavelengths)

Select the spectra to label

n_points = 50

flat_choice = np.random.choice(np.arange(datacube[0].size), n_points, replace=False)
i_points, j_points = np.unravel_index(flat_choice, datacube[0].shape)
np.save('labelled_points.npy', np.array([i_points, j_points]))

i_points, j_points = np.load('labelled_points.npy')

Select the spectra to label from the data file

raw_spectra = datacube[:, i_points, j_points].T

Normalise each spectrum to be in range [0, 1]

labelled_spectra = np.empty((len(raw_spectra), len(model.constant_wavelengths)))
for i in range(len(labelled_spectra)):
    labelled_spectra[i] = normalise_spectrum(raw_spectra[i], model=model)

Script to semi-automate the classification process

• Type a number 0 - 4 for assign a classification to the plotted spectrum

• Type 5 to skip and move on to the next spectrum

• Type back to move to the previous spectrum

• Type exit to give up (keeping ones already done)

The labels are present in the labels variable (-1 represents an unclassified spectrum)

labels = np.full(len(labelled_spectra), -1, dtype=int)
i = 0
while i < len(labelled_spectra):

    # Show the spectrum to be classified along with description
    plt.figure(figsize=(15, 10))
    plt.plot(labelled_spectra[i])
    plt.show()
    print("i = {}".format(i))
    print("absorption --- both --- emission / skip")
    print("       0    1    2    3    4         5 ")

    # Ask for user's classification
    classification = input('Type [0-4]:')

    try:  # Must be an integer
        classification_int = int(classification)
    except ValueError:
        classification_int = -1  # Try current spectrum again

    if classification == 'back':
        i -= 1  # Go back to the previous spectrum
    elif classification == 'exit':
        break  # Exit the loop, saving labels that were given
    elif 0 <= classification_int <= 4:  # Valid classification
        labels[i] = int(classification)  # Assign the classification to the spectrum
        i += 1  # Move on to the next spectrum
    elif classification_int == 5:
        i += 1  # Skip and move on to the next spectrum
    else:  # Invalid integer classification
        i += 0  # Try current spectrum again

Plot bar chart of classification populations

unique, counts = np.unique(labels, return_counts=True)
plt.figure()
plt.bar(unique, counts)
plt.title('Number of spectra in each classification')
plt.xlabel('Classification')
plt.ylabel('N_spectra')
plt.show()

Overplot the spectra of each classification

for classification in unique:
    plt.figure()
    for spectrum in labelled_spectra[labels == classification]:
        plt.plot(model.constant_wavelengths, spectrum)
    plt.title('Classification {}'.format(classification))
    plt.yticks([0, 1])
    plt.show()

Save the labelled spectra for use later

np.save('labelled_data.npy', labelled_spectra)
np.save('labels.npy', labels)

Contributing

Code of Conduct

If you find this package useful and have time to make it even better, you are very welcome to contribute to this package, regardless of how much prior experience you have. Types of ways you can contribute include, expanding the documentation with more use cases and examples, reporting bugs through the GitHub issue tracker, reviewing pull requests and the existing code, fixing bugs and implementing new features in the code. You are encouraged to submit any bug reports and pull requests directly to the GitHub repository.

Please note that this project is released with a Contributor Code of Conduct. By participating in this project you agree to abide by its terms.

Citation

If you have used this package in work that leads to a publication, we would be very grateful if you could acknowledge your use of this package in the main text of the publication. Please cite the following publications,

MacBride CD, Jess DB. 2021 MCALF: Multi-Component Atmospheric Line Fitting. Journal of Open Source Software. 6(61), 3265. (doi:10.21105/joss.03265)

MacBride CD, Jess DB, Grant SDT, Khomenko E, Keys PH, Stangalini M. 2020 Accurately constraining velocity information from spectral imaging observations using machine learning techniques. Philosophical Transactions of the Royal Society A. 379, 2190. (doi:10.1098/rsta.2020.0171)

Please also cite the Zenodo DOI for the package version you used. Please also consider integrating your code and examples into the package.

License

MCALF is licensed under the terms of the BSD 2-Clause license.

Code Reference

MCALF

mcalf Package

MCALF: Multi-Component Atmospheric Line Fitting

MCALF is an open-source Python package for accurately constraining velocity information from spectral imaging observations using machine learning techniques.

MCALF models

This sub-package contains:

• Base and sample models that can be adapted and fitted to any spectral imaging dataset.

• Models optimised for particular data sets that can be used directly.

• Data structures for storing and exporting the fitted parameters, as well as simplifying the calculation of velocities.

mcalf.models Package

Classes

FitResult(fitted_parameters, fit_info)	Class that holds the result of a fit.
FitResults(shape, n_parameters[, time])	Class that holds multiple fit results in a way that can be easily processed.
IBIS8542Model(*[, constant_wavelengths, ...])	Class for working with IBIS 8542 Å calcium II spectral imaging observations.
ModelBase(*[, original_wavelengths, ...])	Base class for spectral line model fitting.

Class Inheritance Diagram

MCALF profiles

This sub-package contains:

• Functions that can be used to model the spectra.

• Voigt profile with a variety of wrappers for different applications (mcalf.profiles.voigt).

• Gaussian profiles and skew normal distributions (mcalf.profiles.gaussian).

mcalf.profiles Package

mcalf.profiles.voigt Module

Functions

voigt_integrate(x, s, g[, clib])	Voigt function implementation (calculated by integrating).
voigt_faddeeva(x, s, g, **kwargs)	Voigt function implementation (Faddeeva).
voigt_mclean(x, s, g, **kwargs)	Voigt function implementation (efficient approximation).
voigt_nobg(x, a, b, s, g[, impl])	Voigt function with no background (Base Voigt function).
voigt(x, a, b, s, g, d, **kwargs)	Voigt function with background.
double_voigt_nobg(x, a1, b1, s1, g1, a2, b2, ...)	Double Voigt function with no background.
double_voigt(x, a1, b1, s1, g1, a2, b2, s2, ...)	Double Voigt function with background.

mcalf.profiles.gaussian Module

Functions

single_gaussian(x, a, b, c, d)

Gaussian function.

MCALF visualisation

This sub-package contains:

• Functions to plot the input spectrum and the fitted model.

• Functions to plot the spatial distribution and their general profile.

• Functions to plot the velocities calculated for a spectral imaging scan.

mcalf.visualisation Package

Functions

bar([class_map, vmin, vmax, reduce, style, ...])	Plot a bar chart of the classification abundances.
init_class_data(class_map[, vmin, vmax, ...])	Initialise dictionary of common classification plotting data.
plot_class_map([class_map, vmin, vmax, ...])	Plot a map of the classifications.
plot_classifications(spectra, labels[, ...])	Plot spectra grouped by their labelled classification.
plot_ibis8542(wavelengths, spectrum[, fit, ...])	Plot an IBIS8542Model fit.
plot_map(arr[, mask, umbra_mask, ...])	Plot a velocity map array.
plot_spectrum(wavelengths, spectrum[, ...])	Plot a spectrum with the wavelength grid shown.

MCALF utils

This sub-package contains:

• Functions for processing spectra (mcalf.utils.spec).

• Functions for smoothing n-dimensional arrays (mcalf.utils.smooth).

• Functions for masking the input data to limit the region computed (mcalf.utils.mask).

• Functions for helping with plotting (mcalf.utils.plot).

• Classes for managing collections of data (mcalf.utils.collections).

• Miscellaneous utility functions (mcalf.utils.misc).

mcalf.utils Package

mcalf.utils.spec Module

Functions

reinterpolate_spectrum(spectrum, ...)	Reinterpolate the spectrum.
normalise_spectrum(spectrum[, ...])	Normalise an individual spectrum to have intensities in range [0, 1].
generate_sigma(sigma_type, wavelengths, ...)	Generate the default sigma profiles.

mcalf.utils.smooth Module

Functions

moving_average(array, width)	Boxcar moving average.
gaussian_kern_3d([width, sigma])	3D Gaussian kernel.
smooth_cube(cube, mask, **kwargs)	Apply Gaussian smoothing to velocities.
mask_classifications(class_map[, vmin, ...])	Mask 2D and 3D arrays of classifications.

mcalf.utils.mask Module

Functions

genmask(width, height[, radius, ...])	Generate a circular mask of specified size.
radial_distances(n_cols, n_rows)	Generates a 2D array of specified shape of radial distances from the centre.

mcalf.utils.plot Module

Functions

hide_existing_labels(plot_settings[, axes, fig])	Hides labels for each dictionary provided if label already exists in legend.
calculate_axis_extent(resolution, px[, ...])	Calculate the extent from a resolution value along a particular axis.
calculate_extent(shape, resolution[, ...])	Calculate the extent from a particular data shape and resolution.
class_cmap(style, n)	Create a listed colormap for a specific number of classifications.

mcalf.utils.collections Module

Classes

Parameter(name[, value])	A named parameter with a optional value.
ParameterDict([dict])	An unordered dictionary of Parameter objects.
OrderedParameterDict	An ordered dictionary of Parameter objects.
BaseParameterDict()	A base class for dictionaries of Parameter objects.
SyncedParameters()	Database for keeping Parameter objects in sync.

Class Inheritance Diagram

mcalf.utils.misc Module

Functions

make_iter(*args)	Returns each inputted argument, wrapping in a list if not already iterable.
load_parameter(parameter[, wl])	Load parameters from file, optionally evaluating variables from strings.
merge_results(filenames, output)	Merges files generated by the mcalf.models.FitResults.save() method.
update_signature(cls)	Update the signature of a model class.

Contributor Covenant Code of Conduct

Our Pledge

We as members, contributors, and leaders pledge to make participation in our community a harassment-free experience for everyone, regardless of age, body size, visible or invisible disability, ethnicity, sex characteristics, gender identity and expression, level of experience, education, socio-economic status, nationality, personal appearance, race, religion, or sexual identity and orientation.

We pledge to act and interact in ways that contribute to an open, welcoming, diverse, inclusive, and healthy community.

Our Standards

Examples of behavior that contributes to a positive environment for our community include:

• Demonstrating empathy and kindness toward other people

• Being respectful of differing opinions, viewpoints, and experiences

• Giving and gracefully accepting constructive feedback

• Accepting responsibility and apologizing to those affected by our mistakes, and learning from the experience

• Focusing on what is best not just for us as individuals, but for the overall community

Examples of unacceptable behavior include:

• The use of sexualized language or imagery, and sexual attention or advances of any kind

• Trolling, insulting or derogatory comments, and personal or political attacks

• Public or private harassment

• Publishing others’ private information, such as a physical or email address, without their explicit permission

• Other conduct which could reasonably be considered inappropriate in a professional setting

Enforcement Responsibilities

Community leaders are responsible for clarifying and enforcing our standards of acceptable behavior and will take appropriate and fair corrective action in response to any behavior that they deem inappropriate, threatening, offensive, or harmful.

Community leaders have the right and responsibility to remove, edit, or reject comments, commits, code, wiki edits, issues, and other contributions that are not aligned to this Code of Conduct, and will communicate reasons for moderation decisions when appropriate.

Scope

This Code of Conduct applies within all community spaces, and also applies when an individual is officially representing the community in public spaces. Examples of representing our community include using an official e-mail address, posting via an official social media account, or acting as an appointed representative at an online or offline event.

Enforcement

Instances of abusive, harassing, or otherwise unacceptable behavior may be reported to the community leaders responsible for enforcement at mcalf@macbride.me. All complaints will be reviewed and investigated promptly and fairly.

All community leaders are obligated to respect the privacy and security of the reporter of any incident.

Enforcement Guidelines

Community leaders will follow these Community Impact Guidelines in determining the consequences for any action they deem in violation of this Code of Conduct:

1. Correction

Community Impact: Use of inappropriate language or other behavior deemed unprofessional or unwelcome in the community.

Consequence: A private, written warning from community leaders, providing clarity around the nature of the violation and an explanation of why the behavior was inappropriate. A public apology may be requested.

2. Warning

Community Impact: A violation through a single incident or series of actions.

Consequence: A warning with consequences for continued behavior. No interaction with the people involved, including unsolicited interaction with those enforcing the Code of Conduct, for a specified period of time. This includes avoiding interactions in community spaces as well as external channels like social media. Violating these terms may lead to a temporary or permanent ban.

3. Temporary Ban

Community Impact: A serious violation of community standards, including sustained inappropriate behavior.

Consequence: A temporary ban from any sort of interaction or public communication with the community for a specified period of time. No public or private interaction with the people involved, including unsolicited interaction with those enforcing the Code of Conduct, is allowed during this period. Violating these terms may lead to a permanent ban.

4. Permanent Ban

Community Impact: Demonstrating a pattern of violation of community standards, including sustained inappropriate behavior, harassment of an individual, or aggression toward or disparagement of classes of individuals.

Consequence: A permanent ban from any sort of public interaction within the community.

Attribution

This Code of Conduct is adapted from the Contributor Covenant, version 2.0, available at https://www.contributor-covenant.org/version/2/0/code_of_conduct.html.

Community Impact Guidelines were inspired by Mozilla’s code of conduct enforcement ladder.

For answers to common questions about this code of conduct, see the FAQ at https://www.contributor-covenant.org/faq. Translations are available at https://www.contributor-covenant.org/translations.

MCALF Licence

MCALF is licensed under the terms of the BSD 2-Clause license.

Copyright (c) 2020 Conor MacBride All rights reserved.

Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:

• Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.

• Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS “AS IS” AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.