{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "9f9d6fc7-af39-42ca-8e85-a83957bf0946",
   "metadata": {},
   "source": [
    "# Coreloss Example\n",
    "In this notebook, a typical workflow of EELS processing is shown on simulated data. This notebook shows one example of how core-loss quantification could be performed "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "0b9c69f6-af2d-4290-a305-9f6c197e7f42",
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib qt \n",
    "#important for the em.MultiSpectrumVisualizer since this is an interactive plotting tool"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "c057e443-1c9a-4a6f-8d76-07403ec6a9c2",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pyEELSMODEL.api as em"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "917c7b16-a202-49f2-bc0c-d910185964b6",
   "metadata": {},
   "source": [
    "## STEM-EELS Simulation\n",
    "A STEM-EELS map is simulated which will be used to quantify the spectra. For more information on how the simulation works one could check the *CoreLossExampleExtended* notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "181bbd08-faf4-4d88-b996-d5a685806239",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from pyEELSMODEL.misc.data_simulator import simulate_data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "22b0e9ad-3682-4221-8324-8028e8de4f49",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "16384it [00:12, 1284.04it/s]\n",
      "2163it [00:04, 457.08it/s]"
     ]
    }
   ],
   "source": [
    "hl, ll = simulate_data()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "38cdfa54-76d4-4088-a169-574f3a72028f",
   "metadata": {},
   "source": [
    "## Quantification Part\n",
    "The previous part is the more irrelevant part which simulates the core-loss. In most of the cases, the data is available and needs to be processed. In this part, a typical workflow on getting the elemental abundance together with estimated errors is shown. Note that the ElementalQuantification class has such a workflow implemented in it but here we show how one can design their own workflow to optimize the data processing."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "475323f7-969f-4e73-9ae7-9e2c31292b17",
   "metadata": {},
   "source": [
    "#### Aligning multispectra"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9f16e0fb-8b8f-4558-9557-61b881d7d5ec",
   "metadata": {},
   "outputs": [],
   "source": [
    "fast_align = em.FastAlignZeroLoss(ll, other_spectra=[hl], cropping=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1fa570a8-0233-4ce2-89ac-1205e7f86cb4",
   "metadata": {},
   "outputs": [],
   "source": [
    "fast_align.perform_alignment()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c80b972c-31b7-4993-89f7-d12e2c913d7d",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = fast_align.show_shift()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b5729c2a-130e-4308-889c-10a4caf15905",
   "metadata": {},
   "outputs": [],
   "source": [
    "hl_al = fast_align.aligned_others[0] #coreloss which is used for quantification\n",
    "ll_al = fast_align.aligned #lowloss which is used for quantification"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b48f4d46-28e2-4b3f-b019-8139349ab858",
   "metadata": {},
   "source": [
    "#### Define the model\n",
    "The next step is to chose a proper model for the experimental data. In pyEELMODEL each model consist out of  components and each component has multiple parameters. For instance, a gaussian component has three parameters: ampltide, center and fwhm. For each parameter one can identify if it needs to be changeable or not. If it is not changeable then it will not be updated via the fitting procedure. The cross sections have multiple parameters but for elemental quantification only the amplitude is the unkown. More information on the model-based approach can be found here [[1]](https://doi.org/10.1016/j.ultramic.2006.05.006) \\\n",
    "The model for this example consists in this case out of three parts:\n",
    "\n",
    "1. **The background**: Historically a powerlaw was used to model the background but in this example a linear background model is used[[2]](https://doi.org/10.1016/j.ultramic.2023.113830). This keeps the entire model linear which is advantages because no starting parameters are needed and no iterations need to be performed to find the global optimum.\n",
    "   \n",
    "2. **Atomic cross sections**: The generalized oscillator strengths from Zhang et al. [[3]](https://zenodo.org/records/7729585) are used. To properly calculate these cross sections, the acceleration voltage (E0), convergence angle (alpha) and collection angle (beta) are needed as input.\n",
    " \n",
    "3. **The low loss**: Due to multiple scattering, the shape of cross sections will be modified and this can be taken into account if the low loss is acquired from the same area. Note that the background will not be convoluted in the model since this is hard to incorporate due to the artifacts arising from the boundaries [[1]](https://doi.org/10.1016/j.ultramic.2006.05.006).\\\n",
    "\\\n",
    "[1] Verbeeck J. et al; Model based quantification of EELS spectra; Ultramicroscopy; 2004; doi:[10.1016/j.ultramic.2006.05.006](https://doi.org/10.1016/j.ultramic.2006.05.006)\\\n",
    "[2] Van den Broek W. et al; Convexity constrains on linear background models for electron energy loss spectra; Ultramicroscopy; 2023; doi:[10.1016/j.ultramic.2023.113830](https://doi.org/10.1016/j.ultramic.2023.113830)\\\n",
    "[3] Zhang Z. et al; Generalised oscillator strength for core-shell excitation by fast electron based on Dirac solutions; Zenodo; 2023; doi:[10.5281/zenodo.7729585](https://zenodo.org/records/7729585)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1f30802e-6459-413d-bbba-ffe14fea6a8f",
   "metadata": {},
   "source": [
    "##### Background component\n",
    "The linear combination of fixed powerlaws where the powers are given by rlist:\\\n",
    "$$bg(E) = \\sum_{i=0}^n A_i E^i$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5eea7d5a-de95-486e-9ea5-dcb6874df70a",
   "metadata": {},
   "outputs": [],
   "source": [
    "from pyEELSMODEL.components.linear_background import LinearBG"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f554ec0c-830d-4444-b1e3-d6476c553d31",
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 4\n",
    "bg = LinearBG(specshape=hl_al.get_spectrumshape(), rlist=np.linspace(1,5,n))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "abd626e3-198a-497f-add0-3249b00a0a0e",
   "metadata": {},
   "source": [
    "##### Cross sections\n",
    "The cross sections are calculated using the cross sections of [Zezhong Zhang](https://zenodo.org/records/7729585). In pyEELSMODEL, the hydrogenic K and L edges and the cross section from [Segger, Guzzinati and Kohl](https://zenodo.org/records/7645765) are also available. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a01805e5-1a89-45fe-b0dc-797d28b1af3e",
   "metadata": {},
   "outputs": [],
   "source": [
    "from pyEELSMODEL.components.CLedge.zezhong_coreloss_edgecombined import ZezhongCoreLossEdgeCombined\n",
    "from pyEELSMODEL.components.CLedge.kohl_coreloss_edgecombined import KohlLossEdgeCombined"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "06770cf4-19bc-4773-be9c-0054511199be",
   "metadata": {},
   "outputs": [],
   "source": [
    "elements = ['C', 'N', 'O', 'Fe']\n",
    "edges = ['K', 'K', 'K', 'L']\n",
    "E0 = 300e3 \n",
    "alpha = 1e-9\n",
    "beta = 20e-3 "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ac9d1146-d7d0-420d-9e43-6f77609f30ab",
   "metadata": {},
   "source": [
    "Showcase the difference between the two different GOS tables. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0ed77fbc-6daf-45c2-b221-820db5f3f014",
   "metadata": {},
   "outputs": [],
   "source": [
    "#can take a bit of time since the cross section is calculated from the tabulated GOS arrays\n",
    "fig, ax = plt.subplots(1,len(elements))\n",
    "for ii in range(len(elements)):\n",
    "    compz = ZezhongCoreLossEdgeCombined(hl_al.get_spectrumshape(), 1, E0, alpha,beta, elements[ii], edges[ii])\n",
    "    compz.calculate() #calculates the cross section with given parameters\n",
    "    \n",
    "    compk = KohlLossEdgeCombined(hl_al.get_spectrumshape(), 1, E0, alpha,beta, elements[ii], edges[ii]) \n",
    "    compk.calculate()\n",
    "\n",
    "    ax[ii].plot(compz.energy_axis, compz.data, label='Zhang')\n",
    "    ax[ii].plot(compk.energy_axis, compk.data, label='Kohl')\n",
    "ax[0].legend()\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "537928c3-e23f-4a39-a600-21a6e3097ae3",
   "metadata": {},
   "outputs": [],
   "source": [
    "#can take a bit of time since the cross section is calculated from the tabulated GOS arrays\n",
    "#can chose which cross section you use \n",
    "comp_elements = []\n",
    "A = 1 #amplitude for cross section, since model is linear this value is not super important. For non-linear fitters the starting value could be important\n",
    "for elem, edge in zip(elements, edges):\n",
    "    comp = ZezhongCoreLossEdgeCombined(hl_al.get_spectrumshape(), 1, E0, alpha,beta, elem, edge)\n",
    "    #comp = KohlLossEdgeCombined(hl_al.get_spectrumshape(), 1, E0, alpha,beta, elem, edge)\n",
    "    comp_elements.append(comp)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "75c780b6-241b-40b7-98ef-edf95d4ad7fd",
   "metadata": {},
   "source": [
    "##### Multiple scattering\n",
    "The calculated components are convolved with the low loss if indicated. For instance, the background component will not be convolved with the lowloss "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7104c1ed-5071-4cc2-9743-6ddf333b5362",
   "metadata": {},
   "outputs": [],
   "source": [
    "from pyEELSMODEL.components.MScatter.mscatterfft import MscatterFFT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d7e064bc-398a-418d-abab-36cf6b146b07",
   "metadata": {},
   "outputs": [],
   "source": [
    "llcomp  = MscatterFFT(hl_al.get_spectrumshape(), ll_al)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "89f70dec-446c-4f48-b73b-9c828b9a05f7",
   "metadata": {},
   "source": [
    "##### Model\n",
    "The model gets created by adding all the different components together into a list. It uses this information to calculate the resulting model and can be used as input for the fitter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f9bbe4b1-44e6-45cd-8395-efcadc268027",
   "metadata": {},
   "outputs": [],
   "source": [
    "components = [bg]+comp_elements+[llcomp]\n",
    "mod = em.Model(hl_al.get_spectrumshape(), components)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b86cb936-c134-47a9-a61c-bba8d2bc35e4",
   "metadata": {},
   "outputs": [],
   "source": [
    "#shows the model with the given paramter values. \n",
    "mod.calculate()\n",
    "mod.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dca7c0a5-2e93-4133-b5b3-77b8ef96ff2c",
   "metadata": {},
   "source": [
    "#### Finding optimal parameters for model\n",
    "Since the model we defined is linear, we can use a weighted linear fitter. The weights are determined from the assumption that the noise is poisson distributed. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fa1180c0-ba37-4914-988a-1e5ba9b591a1",
   "metadata": {},
   "outputs": [],
   "source": [
    "#creates a fit object \n",
    "fit = em.LinearFitter(hl_al, mod, use_weights=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "02ff9f63-799d-4165-bd3b-75b91ffec165",
   "metadata": {},
   "outputs": [],
   "source": [
    "fit.multi_fit()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ebb5a8e4-fc51-48d9-ae8a-fbd029754054",
   "metadata": {},
   "source": [
    "The fitted parameters can be accessed by following functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ec957ff7-11c6-473b-b16e-b9d94588809b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# the fitted parameters can be found in .coeff_matrix attribute.\n",
    "print(fit.coeff_matrix.shape)\n",
    "print(str(fit.coeff_matrix.shape[2])+' parameters are optimized in the fitting procedure')\n",
    "\n",
    "#To know which parameter corresponds to which index in the coeff_matrix, following function can be used\n",
    "N_index = fit.get_param_index(comp_elements[1].parameters[0]) #comp_elements[1].parameters[0]: amplitude of nitrogen edge\n",
    "N_map = fit.coeff_matrix[:,:,N_index]\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.imshow(N_map)\n",
    "ax.set_title(comp_elements[1].name)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b22a467c-baba-4916-bc8e-2830b8b426b2",
   "metadata": {},
   "outputs": [],
   "source": [
    "#function shows the elemental maps \n",
    "fig, maps, names = fit.show_map_result(comp_elements)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ffd5822a-2d62-4d8b-ab13-48a8de73cbb2",
   "metadata": {},
   "outputs": [],
   "source": [
    "#calculates the fitted model, this can be used to validate visually the fitted results\n",
    "multimodel = fit.model_to_multispectrum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ede6cc43-6a12-4b41-ac4f-133f25c6128a",
   "metadata": {},
   "outputs": [],
   "source": [
    "em.MultiSpectrumVisualizer([hl_al, multimodel])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3a7277a6-f30f-4ab1-9dea-8d7e0e1fe6c0",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#calculates the cramer rao lower bound for the given paramter at each probe position\n",
    "crlb_Fe = fit.CRLB_map(comp_elements[3].parameters[0]) #comp_elements[3].parameters[0]: amplitude of iron edge"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "985c54c0-2e71-48b3-89d5-37f462303709",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.imshow(crlb_Fe)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e348463f-e5c0-4d61-9ec6-816e9c3a74de",
   "metadata": {},
   "source": [
    "### ElementalQuantification class\n",
    "The last part shows how the ElementalQuantification class is used as workflow to get the same result. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fc289a1d-8178-49db-bce6-25d9d1277532",
   "metadata": {},
   "outputs": [],
   "source": [
    "settings = (E0, alpha, beta)\n",
    "quant = em.ElementalQuantification(hl, elements, edges, settings, ll=ll)\n",
    "quant.n_bgterms = 4\n",
    "quant.linear_fitter_method = 'ols'\n",
    "quant.do_procedure()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "59ec8335-3942-49f8-ae50-79bd4bf55de6",
   "metadata": {},
   "outputs": [],
   "source": [
    "quant.show_elements_maps()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}