docs/guide/examples/example1/FittingIBIS.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Using the package to fit multi-component atmospheres using the IBIS8542Model class"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Initialising the model and loading the data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import the packages required for the analysis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import mcalf.models\n",
    "import numpy as np\n",
    "from scipy.io import readsav\n",
    "from astropy.io import fits"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load the wavelength positions that IBIS collected data at. Also load the prefilter correction \n",
    "to divide the raw spectra by.\n",
    "Example shows how IDL SAV files can be read along with CSV files."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "original_wavelengths = readsav('wavelengths_original.sav')['wavelengths']\n",
    "prefilter_response = np.loadtxt('prefilter.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Initialise the IBIS 8542 Model with the default parameters and the specific wavelengths and prefilter. \n",
    "Alternatively, a `config` file can be specified."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = mcalf.models.IBIS8542Model(original_wavelengths=original_wavelengths, prefilter_response=prefilter_response)\n",
    "#m = mcalf.models.IBIS8542Model(config=\"config.yml\")  # See `config.yml` for more details."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import the neural network ground truth dataset and train it. Alternativly, a pre-trained neural network can be loaded."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pickle  # Load the trained model from file.\n",
    "# pkl = open('trainedneuralnetwork.pkl', 'rb')\n",
    "# m.neural_network = pickle.load(pkl)  # Overwrite the default untrained model\n",
    "\n",
    "# Load ground truth data\n",
    "labels = np.load('labels.npy')\n",
    "spectra = np.load('labeled_data.npy')\n",
    "\n",
    "# Train the neural network on the first half of the ground truth data\n",
    "m.train(spectra[:100], labels[:100])\n",
    "\n",
    "# Optionally, save the trained neural network so it can be loaded again later\n",
    "# pkl = open('trainedneuralnetwork.pkl', 'wb')\n",
    "# pickle.dump(m.neural_network, pkl)\n",
    "# pkl.close()\n",
    "\n",
    "# Test the neural network with the rest of the ground truth data\n",
    "m.test(spectra[100:], labels[100:])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load the spectral data from file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hdul = fits.open('IBIS_scan_00100.fits')\n",
    "raw = hdul[0].data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Mask the region that you would like to process."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load the mask\n",
    "from mcalf.utils.mask import genmask\n",
    "mask = genmask(width=1000, height=1000, radius=460, right_shift=-5, up_shift=10)\n",
    "# mask = np.load('umbra_mask.npy')  # Or use a premade one that can be any shape.\n",
    "\n",
    "# Apply the mask\n",
    "mask = np.repeat(mask[np.newaxis, :, :], len(raw), axis=0)  # Duplicate along the wavelength dimension of `raw'\n",
    "raw[~mask] = np.nan  # Code will ignore nan pixels"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import the masked data into the IBIS 8542 Model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m.load_array(raw, ['wavelength', 'row', 'column'])\n",
    "\n",
    "# Multiple times can be loaded at once, but you'll have to adapt the above masking code\n",
    "# m.load_array(umbra_data, ['time', 'wavelength', 'row', 'column'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load the calculated background values (continuum intensities) for all the spectra."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bg = np.load('background_averaged.npy', mmap_mode='r')[100] \n",
    "# Use mmap_mode so you can quickly pick the scan that you are working with without loading the whole file\n",
    "m.load_background(bg, ['row', 'column'])\n",
    "\n",
    "# Multiple times can be loaded at once\n",
    "# m.load_background(bg, ['time', 'row', 'column'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now the model has been fully initialised. You may wish to adjust some more parameters. \n",
    "See IBIS8542Model docstring for details."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using the model to classify spectra"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The loaded spectra can be now classified."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Select the data that you would like to classify\n",
    "classifications_map = m.classify_spectra(row=range(1000), column=range(1000))\n",
    "# classifications_map = m.classify_spectra(row=[300, 301, 302], column=range(500, 505))\n",
    "\n",
    "# The following will give you (300, 500), (300, 600), (400, 500), and (400, 600)\n",
    "# classifications_map = m.classify_spectra(row=[300, 400], column=[500, 600])\n",
    "\n",
    "# Process the result to make it easier to plot\n",
    "classifications_map = np.asarray(classifications_map[0], dtype=float)  # Zero index for the first (and only) loaded time\n",
    "classifications_map[classifications_map == -1] = np.nan  # Replace the invalid values with nan to aid visualisation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can then run analysis on the classifications or plot them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.imshow(classifications_map)\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using the model to fit the spectra and find velocities"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can simply fit a single spectrum and plot it as follows,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fit = m.fit(row=600, column=600)[0]  # The zero index selects the first (and only) FitResult generated\n",
    "\n",
    "# You can then call the plot method on the FitResult, remembering to specift the model used such that it can also plot \n",
    "# the observed spectrum and average central wavelength\n",
    "fit.plot(m)\n",
    "fit.plot(m, separate=True)  # If a double Voigt model is fitted, this will show the two components separately\n",
    "\n",
    "# Alternative, equivalent, method\n",
    "m.plot(fit)\n",
    "m.plot_separate(fit)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you have multiple fits, they will always be returned as a 1D list of FitResult objects, \n",
    "but the you can find the actual index as follows,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fit1, fit2 = m.fit(row=[500, 501], column=500)\n",
    "\n",
    "# Index of fit1 [<time>, <row>, <column>], and other details\n",
    "print(\"Index:\", fit1.index, fit2.index)\n",
    "print(\"Fitted parameters:\", fit1.parameters, fit2.parameters)\n",
    "print(\"Classification assigned:\", fit1.classification, fit2.classification)\n",
    "print(\"Profile used:\", fit1.profile, fit2.profile)\n",
    "print(\"Was the method able to produce a result?\", fit1.success, fit2.success)\n",
    "\n",
    "# You can find the velocity using:\n",
    "\n",
    "#   (the model `m' must be specified such that the stationary line core wavelength is known)\n",
    "# Quiescent is the default velocity if not specified:\n",
    "print(\"Quiescent velocity:\", fit1.velocity(m, vtype='quiescent'), fit2.velocity(m, vtype='quiescent'))\n",
    "# The following will work if a double Voigt fit is done (otherwise nan):\n",
    "print(\"Active velocity:\", fit1.velocity(m, vtype='active'), fit2.velocity(m, vtype='active'))\n",
    "\n",
    "fit1.plot(m)\n",
    "fit2.plot(m)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Running big jobs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you are fitting a very large number of velocities the FitResults class is handy as you can append a FitResult to it \n",
    "and it will extract the data into arrays that can easily be saved."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results = mcalf.models.FitResults((1000, 1000), 8)  # Assuming a loaded array of 1000 by 1000 rows and columns, \n",
    "# and max 8 fitted parameters per spectrum."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can then run the following code to add the previous fits to it,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results.append(fit)\n",
    "results.append(fit1)\n",
    "results.append(fit2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And then you can extract the arrays,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# results.parameters\n",
    "# results.classifications\n",
    "# results.profile\n",
    "# results.success\n",
    "\n",
    "# The velocities can be generated in bulk too\n",
    "\n",
    "# results.velocities(m, vtype='quiescent')\n",
    "# results.velocities(m, vtype='active')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When running a big job you can take the initiallised model, m, from the first section, and run the following,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Assuming 1000 x 1000 rows and columns\n",
    "results = mcalf.models.FitResults((1000, 1000), 8)\n",
    "\n",
    "for i in range(0, 1000, 100):  # Do 10 batches of 100 rows\n",
    "    \n",
    "    # Fit the batch, using the specified number of processing pools\n",
    "    results_batch = m.fit(row=range(i, i+100), column=range(1000), n_pools=32)\n",
    "    \n",
    "    for result in results_batch:  # Take each FitResult object in the batch\n",
    "        results.append(result)  # And append it to the FitResults object"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The arrays can be extracted as shown above, and then saved in your desired format. \n",
    "\n",
    "Any spectra that were masked will be skipped and will not be outputted by the m.fit() function.\n",
    "\n",
    "If you are processing multiple scans, these can be read into the model, m, yet you may have to adapt the code \n",
    "to stay within the memory limitations of your machine. Such workarounds include, \n",
    "calling m.load_array() and m.load_background() multiple time in the same program, \n",
    "i.e. load one IBIS scan, fit it, extract the result, load the next IBIS scan and repeat."
   ]
  }
 ],
 "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}