PK |VnMK K models/plot_ibis8542data.ipynb{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n# Working with IBIS data\nThis example shows how to initialise the\n:class:`mcalf.models.IBIS8542Model` class with\nreal IBIS data, and train a neural network\nclassifier. We then proceed to fit the array\nof spectra and visualise the results.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Download sample data\n\nFirst, the sample data needs to be downloaded\nfrom the GitHub repository where it is hosted.\nThis will create four new files in the current\ndirectory (about 651 KB total).\n**You may need to install the requests\nPython package for this step to run.**\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import requests\n\npath = 'https://raw.githubusercontent.com/ConorMacBride/mcalf/main/examples/data/ibis8542data/'\n\nfor file in ('wavelengths.txt', 'spectra.fits',\n 'training_data.json', 'results.fits'):\n r = requests.get(path + file, allow_redirects=True)\n with open(file, 'wb') as f:\n f.write(r.content)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Load the sample data\n\nNext, the downloaded data needs to be loaded into Python.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Import the packages needed for loading data\nimport json\nimport numpy as np\nfrom astropy.io import fits\n\n# Load the spectra's wavelength points\nwavelengths = np.loadtxt('wavelengths.txt', dtype='>f4')\n\n# Load the array of spectra\nwith fits.open('spectra.fits') as hdul:\n spectra = hdul[0].data\n\n# Load indices of labelled spectra\nwith open('training_data.json', 'r') as f:\n data = f.read()\ntraining_data = json.loads(data)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As you can see, the sample data consists of\na 60 by 50 array of spectra with 27 wavelength\npoints,\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "print(wavelengths.shape, spectra.shape)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The blue wing and line core intensity values\nof the spectra are plotted below for illustrative\npurposes,\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import matplotlib.pyplot as plt\nimport astropy.units as u\nfrom mcalf.visualisation import plot_map\n\nfig, ax = plt.subplots(1, 2, sharey=True, constrained_layout=True)\n\nwing_data = np.log(spectra[0])\ncore_data = np.log(spectra[len(wavelengths)//2])\n\nres = {\n 'offset': (-25, -30),\n 'resolution': (0.098 * 5 * u.arcsec, 0.098 * 5 * u.arcsec),\n 'show_colorbar': False,\n}\n\nwing = plot_map(wing_data, ax=ax[0], **res,\n vmin=np.min(wing_data), vmax=np.max(wing_data))\ncore = plot_map(core_data, ax=ax[1], **res,\n vmin=np.min(core_data), vmax=np.max(core_data))\n\nwing.set_cmap('gray')\ncore.set_cmap('gray')\n\nax[0].set_title('blue wing')\nax[1].set_title('line core')\nax[1].set_ylabel('')\n\nplt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Generate the backgrounds array\n\nAs discussed in the `Using IBIS8542Model`\nexample, a background intensity value\nmust be specified for each spectrum.\n\nFor this small sample dataset, we shall\nsimply use the average of the four leftmost\nintensity values of each spectrum,\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "backgrounds = np.mean(spectra[:4], axis=0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Initialise the IBIS8542Model\n\nThe loaded data can now be passed into\nan :class:`mcalf.models.IBIS8542Model`\nobject.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import mcalf.models\n\nmodel = mcalf.models.IBIS8542Model(original_wavelengths=wavelengths, random_state=0)\n\nmodel.load_background(backgrounds, ['row', 'column'])\nmodel.load_array(spectra, ['wavelength', 'row', 'column'])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Training the neural network\n\nBy default, the :class:`mcalf.models.IBIS8542Model`\nobject is loaded with an untrained neural network,\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "model.neural_network" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The :class:`mcalf.models.IBIS8542Model`\nclass provides two methods to train and test\nthe loaded neural network.\n\nThe ``training_data.json`` file contains a\ndictionary of indices for each classification\n0 to 5. These indices correspond to randomly\npre-chosen spectra in the ``spectra.fits`` file.\n\nThe training set consists of 200 spectra;\n40 for each classification. This training set\nis for demonstration purposes only, generally\nit is not recommended to train with such a\nrelatively high percentage of your data,\nas the risk of overfitting the neural network\nto this specific 60 by 50 array is increased.\n\nTo begin, we'll convert the list of indices\ninto a list of spectra and corresponding\nclassifications,\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from mcalf.utils.spec import normalise_spectrum\n\n\ndef select_training_set(indices, model):\n for c in sorted([int(i) for i in indices.keys()]):\n i = indices[str(c)]\n spectra = np.array([normalise_spectrum(\n model.get_spectra(row=j, column=k)[0, 0, 0],\n model.constant_wavelengths, model.constant_wavelengths\n ) for j, k in i])\n try:\n _X = np.vstack((_X, spectra))\n _y = np.hstack((_y, [c] * len(spectra)))\n except NameError:\n _X = spectra\n _y = [c] * len(spectra)\n return _X, _y\n\n\nX, y = select_training_set(training_data, model)\n\nprint(X.shape) # spectra\nprint(y.shape) # labels/classifications" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "These classifications look as follows,\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from mcalf.visualisation import plot_classifications\n\nplot_classifications(X, y)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we can train the neural network on\n100 labelled spectra (even indices),\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "model.train(X[::2], y[::2])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And now we can use the other 100\nlabelled spectra (odd indices)\nto test the performance of the neural network,\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "model.test(X[1::2], y[1::2])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now that we have a trained neural network,\nwe can use it to classify spectra.\nUsually spectra will be classified automatically\nduring the fitting process, however, you\ncan request the classification by themselves,\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "classifications = model.classify_spectra(row=range(60), column=range(50))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "These classifications look as follows,\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from mcalf.visualisation import plot_class_map\n\nplot_class_map(classifications)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Creating a reproducible classifier\n\nThe neural network classifier introduces a certain amount\nof randomness when it it fitting based on the training\ndata. This randomness arises in the initial values\nof the weights and biases that are fitted during the\ntraining process, as well as the order in which the\ntraining data are used.\n\nThis means that two neural networks trained on identical\ndata will not produce the same results. To aid the\nreproducibility of results that rely on a neural\nnetwork's classifications, a `random_state` integer\ncan be passed to :class:`mcalf.models.IBIS8542Model`\nas we did above. When we set this value to an integer,\nno matter how many times we train the neural network\non the same data, it will always give the same\nresults.\n\nUntil better solutions are available to store trained\nneural networks, a trained neural network can be saved\nto a Python pickle file and later reloaded. For\nmaximum compatibility, it is recommended to reload\ninto the same version of scikit-learn and its\ndependencies.\n\nThe neural network trained above can be saved as follows,\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import pickle\npkl = open('trained_neural_network.pkl', 'wb')\npickle.dump(model.neural_network, pkl)\npkl.close()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This trained neural network can then be reloaded at a\nlater date as follows,\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import pickle\npkl = open('trained_neural_network.pkl', 'rb')\nmodel.neural_network = pickle.load(pkl) # Overwrite the default untrained model" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And you can see that the classifications of spectra are the same,\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "plot_class_map(model.classify_spectra(row=range(60), column=range(50)))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Please see the\n[scikit-learn documentation](https://scikit-learn.org/stable/modules/model_persistence.html)\nfor more details on model persistence.\n\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Fitting the spectra\n\nNow that the data have been loaded and the\nneural network has been trained, we can proceed\nto fit the spectra.\n\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Using pre-calculated results\n\nAs our 60 by 50 array contains 3000 spectra, it would\ntake roughly 10 minutes to fit them all over 6 processing\npools. We include pre-calculated results in the\ndownloaded ``results.fits`` file.\n\nThe next step of the example loads this file back into\nPython as though we have just directly calculated it.\nThis isn't something you would usually need to do,\nso do not worry about the contents of the\n``load_results()`` function, however, we plan to\ninclude this functionality in MCALF itself in the future.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def load_results(file):\n\n with fits.open(file) as hdul:\n for hdu in hdul:\n if hdu.name == 'PARAMETERS':\n r_parameters = hdu.data.copy().reshape(-1, 8)\n elif hdu.name == 'CLASSIFICATIONS':\n r_classifications = hdu.data.copy().flatten()\n r_profile = np.full_like(r_classifications, 'both', dtype=object)\n r_profile[r_classifications <= 1] = 'absorption'\n elif hdu.name == 'SUCCESS':\n r_success = hdu.data.copy().flatten()\n elif hdu.name == 'CHI2':\n r_chi2 = hdu.data.copy().flatten()\n\n results = []\n for i in range(len(r_parameters)):\n fitted_parameters = r_parameters[i]\n fit_info = {\n 'classification': r_classifications[i],\n 'profile': r_profile[i],\n 'success': r_success[i],\n 'chi2': r_chi2[i],\n 'index': [0, *np.unravel_index(i, (60, 50))],\n }\n if fit_info['profile'] == 'absorption':\n fitted_parameters = fitted_parameters[:4]\n results.append(mcalf.models.FitResult(fitted_parameters, fit_info))\n\n return results\n\n\nresult_list = load_results('results.fits')\n\nresult_list[:4] # The first four" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Using your own results\n\nYou can run the following code to generate\nthe ``result_list`` variable for yourself.\nTry starting with a smaller range of rows\nand columns, and set the number of pools\nbased on the specification of your\nprocessor.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# result_list = model.fit(row=range(60), column=range(50), n_pools=6)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The order of the :class:`mcalf.models.FitResult`\nobjects in this list will also differ as\nthe order that spectra finish fitting in\neach pool is unpredictable.\n\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Merging the FitResult objects\n\nThe list of :class:`mcalf.models.FitResult`\nobjects can be mergerd into a\n:class:`mcalf.models.FitResults` object,\nand then saved to a file, just like the\n``results.fits`` file downloaded earlier.\n\nFirst the object needs to be initialised\nwith the spatial dimensions and number\nof fitted parameters,\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "results = mcalf.models.FitResults((60, 50), 8)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we can loop through the list of fits\nand append them to this object,\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "for fit in result_list:\n results.append(fit)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This object can then be saved to file,\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "results.save('ibis8542data.fits', model)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The file has the following structure,\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "with fits.open('ibis8542data.fits') as hdul:\n hdul.info()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exploring the fitted data\n\nYou can plot a fitted spectrum as follows,\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "model.plot(result_list[0])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "You can calculate and plot Doppler velocities\nfor both the quiescent and active regimes\nas follows, (with an outline of the sunspot's\numbra),\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "vq = results.velocities(model, vtype='quiescent')\nva = results.velocities(model, vtype='active')\n\numbra_mask = np.full_like(backgrounds, True)\numbra_mask[backgrounds < 1100] = False\n\nfig, ax = plt.subplots(1, 2, sharey=True, constrained_layout=True)\n\nsettings = {\n 'show_colorbar': False, 'vmax': 4, 'offset': (-25, -30),\n 'resolution': (0.098 * 5 * u.arcsec, 0.098 * 5 * u.arcsec),\n 'umbra_mask': umbra_mask,\n}\n\nim = plot_map(vq, ax=ax[0], **settings)\nplot_map(va, ax=ax[1], **settings)\n\nax[0].set_title('quiescent')\nax[1].set_title('active')\nax[1].set_ylabel('')\nfig.colorbar(im, ax=ax, location='bottom', label='Doppler velocity (km/s)')\n\nplt.show()" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.6" } }, "nbformat": 4, "nbformat_minor": 0 }PK |V"7; 7; models/plot_ibis8542model.ipynb{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n# Using IBIS8542Model\nThis is an example showing how to fit an array of\nspectra using the :class:`mcalf.models.IBIS8542Model`\nclass, and plot and export the results.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Generate sample data\n\n### Create spectral grid\n\nFirst, we shall generate some random sample data to\ndemonstrate the API.\nThe randomly generated spectra are not intended to\nbe representative of observed spectra numerically,\nthey just have a similar shape.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import numpy as np\nimport matplotlib.pyplot as plt\n\nfrom mcalf.profiles.voigt import voigt, double_voigt\nfrom mcalf.visualisation import plot_spectrum\n\n# Create demo wavelength grid\nw = np.linspace(8541, 8543, 20)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We shall generate demo spectra in a 2x3 grid.\nHalf of the spectra will have one spectral\ncomponent (absorption only) and the other\nhalf will have two spectral components\n(mixed absorption and emission components).\n\nDemo spectral components will be modelled as\nVoigt functions with randomly generated\nparameters, each within a set range of values.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def v(classification, w, *args):\n \"\"\"Voigt function wrapper.\"\"\"\n if classification < 1.5: # absorption only\n return voigt(w, *args[:4], args[-1])\n return double_voigt(w, *args) # absorption + emission\n\n\ndef s():\n \"\"\"Generate random spectra for a 2x3 grid.\"\"\"\n np.random.seed(0) # same spectra every time\n p = np.random.rand(9, 6) # random Voigt parameters\n # 0 1 2 3 4 5 6 7 8 # p index\n # a1 b1 s1 g1 a2 b2 s2 g2 d # Voigt parameter\n # absorption |emission |background\n\n p[0] = 100 * p[0] - 1000 # absorption amplitude\n p[4] = 100 * p[4] + 1000 # emission amplitude\n\n for i in (1, 5): # abs. and emi. peak positions\n p[i] = 0.05 * p[i] - 0.025 + 8542\n\n for i in (2, 3, 6, 7): # Voigt sigmas and gammas\n p[i] = 0.1 * p[i] + 0.1\n\n p[8] = 300 * p[8] + 2000 # intensity background constant\n\n # Define each spectrum's classification\n c = [0, 2, 0, 2, 0, 2]\n # Choose single or double component spectrum\n # based on this inside the function `v()`.\n\n # Generate the spectra\n specs = [v(c[i], w, *p[:, i]) for i in range(6)]\n\n # Reshape to 2x3 grid\n return np.asarray(specs).reshape((2, 3, len(w)))\n\n\nraw_data = s()\n\nprint('shape of spectral grid:', raw_data.shape)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "These spectra look as follows,\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fig, axes = plt.subplots(2, 3, constrained_layout=True)\nfor ax, spec in zip(axes.flat, raw_data.reshape((6, raw_data.shape[-1]))):\n\n plot_spectrum(w, spec, normalised=False, ax=ax)\n\nplt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Compute the background intensity\n\n`MCALF` does not model a constant background value,\ni.e., the fitting process assumes the intensity\nvalues far out in the spectrum's wings are zero.\n\nYou can however tell `MCALF` what the constant\nbackground value is, and it will subtract it\nfrom the spectrum before fitting.\n\nThis process was not made automatic as we wanted\nto give the user full control on setting the\nbackground value.\n\nFor these demo data, we shall simply set the background\nto the first intensity value of every spectrum.\nFor a real dataset, you may wish to average the\nbackground value throughout a range of spectral points,\nor even do a moving average throughout time.\nFunctions are provided in ``mcalf.utils.smooth`` to\nassist with this process.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "backgrounds = raw_data[:, :, 0]\n\nprint('shape of background intensity grid:', backgrounds.shape)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Define a demo classifier\n\nIn this example we will not demonstrate how to create\na neural network classifier. `MCALF` offers a lot of\nflexibility when it comes to the classifier. This\ndemo classifier outlines the basic API that is\nrequired. By default, the model has a\n:class:`sklearn.neural_network.MLPClassifier`\nobject preloaded for use as the classifier.\nThe methods :meth:`mcalf.models.ModelBase.train`\nand :meth:`mcalf.models.ModelBase.test` are\nprovided by `MCALF` to assist with training the\nneural network. There is also a useful script\nin the `Getting Started` section under\n`User Documentation` in the sidebar for\nsemi-automating the process of creating the\nground truth dataset.\n\nPlease see tutorials for packages such\nas scikit-learn for more in-depth advice on\ncreating classifiers.\n\nAs we only have six spectra with two distinct\nshapes, we can create a very simple classifier\nthat classifies spectra based on whether their\ncentral intensity is greater or smaller than\nthe left most intensity.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "class DemoClassifier:\n @staticmethod\n def predict(X):\n y = np.zeros(len(X), dtype=int)\n y[X[:, len(X[0]) // 2] > X[:, 0]] = 2\n return y" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Using ``IBIS8542Model`` with the generated data\n\n### Initialise the model\n\nEverything we have been doing up to this point has been\ncreating the demo data and classifier. Now we can actually\ncreate a model and load in the demo data, although the\nfollowing steps would be identical for a real dataset.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from mcalf.models import IBIS8542Model\n\n# Initialise the model with the wavelength grid\nmodel = IBIS8542Model(original_wavelengths=w)\n\n# Load the spectral shape classifier\nmodel.neural_network = DemoClassifier\n\n# Load the array of spectra and background intensities\nmodel.load_array(raw_data, ['row', 'column', 'wavelength'])\nmodel.load_background(backgrounds, ['row', 'column'])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Fit the loaded data\n\nWe have now fully initialised the model. We can now\ncall methods to fit the model to the loaded spectra.\nIn the following step we fit all the loaded spectra\nand a 1D list of :class:`~mcalf.models.FitResult`\nobjects is returned.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fits = model.fit(row=[0, 1], column=range(3))\n\nprint(fits)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Plot the fitted parameters\n\nThe individual components of each fit can now be\nplotted separately,\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fig, axes = plt.subplots(2, 3, constrained_layout=True)\nfor ax, fit in zip(axes.flat, fits):\n\n model.plot_separate(fit, show_legend=False, ax=ax)\n\nplt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Extract processed spectra\n\nAs well as fitting spectra, we can call other methods.\nIn this step we'll extract the array of loaded spectra.\nHowever, these spectra have been re-interpolated to a\nnew finer wavelength grid. (This grid can be customised\nwhen initialising the model.)\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "spectra = model.get_spectra(row=[0, 1], column=range(3))\n\nprint('new shape of spectral grid:', spectra.shape)\n\nspectra_1d = spectra.reshape((6, spectra.shape[-1]))\n\nfig, axes = plt.subplots(2, 3, constrained_layout=True)\nfor ax, spec in zip(axes.flat, spectra_1d):\n\n plot_spectrum(model.constant_wavelengths, spec,\n normalised=False, ax=ax)\n\nplt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Classify spectra\n\nWe can also classify the loaded spectra and create plots,\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "classifications = model.classify_spectra(row=[0, 1], column=range(3))\n\nprint('classifications are:', classifications)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This function plots the spectra grouped\nby classification,\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from mcalf.visualisation import plot_classifications\n\nplot_classifications(spectra_1d, classifications.flatten())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This function plots a spatial map of\nthe classifications on the 2x3 grid,\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from mcalf.visualisation import plot_class_map\n\nplot_class_map(classifications)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Merge output data\n\nThe :class:`~mcalf.models.FitResult` objects in the\n``fits`` 1D list can be merged into a grid.\nEach of these objects can be appended to a\nsingle :class:`~mcalf.models.FitResults` object.\nThis allows for increased data portability,\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from mcalf.models import FitResults\n\n# Initialise with the grid shape and num. params.\nresults = FitResults((2, 3), 8)\n\n# Append each fit to the object\nfor fit in fits:\n results.append(fit)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now that we have appended all 6 fits,\nwe can access the merged data,\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "print(results.classifications)\nprint(results.profile)\nprint(results.success)\nprint(results.parameters)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Calculate velocities\n\nAnd finally, we can calculate Doppler\nvelocities for both the quiescent (absorption)\nand active (emission) regimes.\n(The model needs to be given so the stationary\nline core wavelength is available.)\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "quiescent = results.velocities(model, vtype='quiescent')\nactive = results.velocities(model, vtype='active')\n\nfrom mcalf.visualisation import plot_map\n\nfig, axes = plt.subplots(2, constrained_layout=True)\nplot_map(quiescent, ax=axes[0])\nplot_map(active, ax=axes[1])\nplt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Export to FITS file\n\nThe :class:`~mcalf.models.FitResults` object can\nalso be exported to a FITS file,\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "results.save('ibis8542model_demo_fit.fits', model)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The file has the following structure,\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from astropy.io import fits\n\nwith fits.open('ibis8542model_demo_fit.fits') as hdul:\n hdul.info()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This example outlines the basics of using the\n:class:`mcalf.models.IBIS8542Model` class.\nTypically millions of spectra would need to be\nfitted and processed at the same time.\nFor more details on running big jobs and how\nspectra can be fitted in parallel, see\nthe `Working with IBIS data` example in the\ngallery.\n\n
Due to limitations with the documentation\n hosting provider, the examples in the MCALF\n documentation are computed using a Python\n based implementation of the Voigt profile\n instead of the more efficient C version.\n When this code is run on your own machine\n it should take much less time than the\n value reported below.