{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "d3a83569-6e57-431e-ac91-227a4b5d43a6",
   "metadata": {},
   "source": [
    "# Notebook showing the ML-MGM method implemented in pyFRESCO"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "64fa6b08-95d6-4089-914a-62a29f99dfde",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pyfresco"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4b61083f-0ce3-4983-b702-8c228640b33f",
   "metadata": {},
   "source": [
    "## Data import\n",
    "\n",
    "In the following code cell we define the paths pointing to the img and header files of the proper hyperspectral datacube and of the summary parameters datacube.\n",
    "These two datacubes are the ones needed by pyFRESCO to work properly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "fade6b0e-6d27-41ad-a2a7-d3e547802a00",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['R770', 'RBR', 'BD530_2', 'SH600_2', 'SH770', 'BD640_2', 'BD860_2', 'BD920_2', 'RPEAK1', 'BDI1000VIS', 'R440', 'IRR1', 'BDI1000IR', 'OLINDEX3', 'R1330', 'BD1300', 'LCPINDEX2', 'HCPINDEX2', 'VAR', 'ISLOPE1', 'BD1400', 'BD1435', 'BD1500_2', 'ICER1_2', 'BD1750_2', 'BD1900_2', 'BD1900R2', 'BDI2000', 'BD2100_2', 'BD2165', 'BD2190', 'MIN2200', 'BD2210_2', 'D2200', 'BD2230', 'BD2250', 'MIN2250', 'BD2265', 'BD2290', 'D2300', 'BD2355', 'SINDEX2', 'ICER2_2', 'MIN2295_2480', 'MIN2345_2537', 'BD2500_2', 'BD3000', 'BD3100', 'BD3200', 'BD3400_2', 'CINDEX2', 'BD2600', 'IRR2', 'IRR3', 'R530', 'R600', 'R1080', 'R1506', 'R2529', 'R3920']\n"
     ]
    }
   ],
   "source": [
    "path_if_mtrdr = \"FRT00009b5a/frt00009b5a_07_if165j_mtr3.img\" # Path to the proper hyperspectral datacube\n",
    "head_if_mtrdr = \"FRT00009b5a/frt00009b5a_07_if165j_mtr3.hdr\" # Path to the header of the proper hyperspectral datacube\n",
    "\n",
    "path_sr_mtrdr = \"FRT00009b5a/frt00009b5a_07_sr165j_mtr3.img\" # Path to the summary parameters datacube\n",
    "head_sr_mtrdr = \"FRT00009b5a/frt00009b5a_07_sr165j_mtr3.hdr\" # Path to the header of the summary parameters datacube\n",
    "\n",
    "img , img_sr , wavelength , sr_names = pyfresco.open_raw(path_if_mtrdr , head_if_mtrdr , path_sr_mtrdr , head_sr_mtrdr) # Open the two datacubes\n",
    "\n",
    "Nbands = len(wavelength)\n",
    "print(sr_names)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "476af2cd-a7eb-4a9d-9a6e-a6351c8367a5",
   "metadata": {},
   "source": [
    "## Input spectrum\n",
    "\n",
    "The ML-MGM method is applied to a normalized spectrum.\n",
    "Here we reproduce a compact extraction and normalization workflow to create a target spectrum."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "0ca73990-0e80-49a4-9d09-3015eab0f5a0",
   "metadata": {},
   "outputs": [],
   "source": [
    "SE = pyfresco.SpectraExtract(img , Nbands , wavelength , 750 , 2600) # Definition of the object for spectral extraction\n",
    "\n",
    "RGB = SE.upload_map('MAF' , folder = 'FRT00009b5a/') # Upload an RGB map for the ROI selection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "be23cffe-f04b-4c91-8b54-554080a87c93",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of points inside the drewn polygon:  425\n",
      "64.25992779783394 % of the errors, set to zero for simplicity, are in reality NaN values.\n"
     ]
    }
   ],
   "source": [
    "target_spectra , L_target , target_mask = SE.polygon_spectra(save_pixel = True ,\n",
    " folder = 'FRT00009b5a' ,\n",
    " name = 'mgm_target_coordinates') # Select the target ROI\n",
    "\n",
    "target_spectrum , target_error = SE.final_spectra(mean = True , c = 'blue') # Target mean spectrum\n",
    "\n",
    "SN = pyfresco.SpectraNorm(RGB , Nbands , img , img_sr , wavelength ,\n",
    "                         target_spectrum , target_error , target_spectra ,\n",
    "                         750 , 2600)\n",
    "\n",
    "neutral_spectra , neutral_spectrum , hull_error = SN.neutral_convex_hull() # Select the neutral ROI\n",
    "\n",
    "norm_mgm , error_mgm = SN.norm_spectra(convex_hull=True) # Normalized spectrum\n",
    "SN.normplot()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8fa79abb-2132-43c9-8c1b-ca6629cc5980",
   "metadata": {},
   "source": [
    "The output of the above cell should look like this: \n",
    "![targetmgmroi](notebook5_images/ROI_selection.png)\n",
    "![convexhullroi](notebook5_images/convex_hull.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "225cd277-8b6e-4f3f-9c31-d16a58a6629b",
   "metadata": {},
   "source": [
    "## Definition of the SpectraAnalysis object\n",
    "\n",
    "The MGM routine is part of the `SpectraAnalysis` class.\n",
    "The object is initialized with the normalized spectrum and with the spectra used to obtain it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "id": "c475eead-2055-4273-8d6d-e29e70f3c50c",
   "metadata": {},
   "outputs": [],
   "source": [
    "SA = pyfresco.SpectraAnalysis(norm_mgm , error_mgm ,\n",
    " m_spec = target_spectrum , err_spec = target_error ,\n",
    " n_spectra = neutral_spectrum , n_err = neutral_error ,\n",
    " wavelength = wavelength ,\n",
    " MIN = 750 , MAX = 2600 ,\n",
    " folder = None) # Definition of the object for spectral analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8029ded8-53b7-4439-8543-2c51fe5113aa",
   "metadata": {},
   "source": [
    "## ML-MGM with fixed band centers\n",
    "\n",
    "If the approximate absorption positions or their bounds are known, they can be called as parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "id": "ef7cdc84-39a0-40a0-beb5-d8b27e96d311",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Initial Conditions are: \n",
      " A =  [0.2791449  0.24437045 0.13631059] \n",
      " FWHM =  [ 750.13153  994.27216 1288.158  ] \n",
      " Alpha =  [-1.0107338 -0.9217194 -1.0085499] \n",
      " M =  [1704.6702 1573.8821 1158.6619] \n",
      "\n"
     ]
    }
   ],
   "source": [
    "losses , gaussian_sums , n_wav , gaussian_sum , a , m , fwhm , alpha = SA.mgm(\n",
    "                                                                         n = 3 ,\n",
    "                                                                         iterations = 25000 ,\n",
    "                                                                         lr = 0.001 ,\n",
    "                                                                         means = None,\n",
    "                                                                         mean_ranges = None,\n",
    "                                                                         asym = True ,\n",
    "                                                                         smooth = True ,\n",
    "                                                                         hull = True ,\n",
    "                                                                         t = 'savgol' ,\n",
    "                                                                         w = 7 ,\n",
    "                                                                         o = 2) # Fit four skew-normal absorptions"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6650e904-ff26-4881-9c40-5999a0e123f9",
   "metadata": {},
   "source": [
    "## Fit diagnostics\n",
    "\n",
    "The following methods evaluate and visualize the global fit, the local contribution of each band and the residuals."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "id": "05f08085-56f9-4b79-9aea-c660b6cb4053",
   "metadata": {},
   "outputs": [],
   "source": [
    "residuals = SA.mse() # Mean square error per wavelength\n",
    "\n",
    "SA.global_fit() # Global fit and residuals\n",
    "\n",
    "single_bands = SA.local_fit() # Single fitted distributions\n",
    "\n",
    "all_bands = SA.both_fit() # Global fit, local bands and residuals"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b1dc4b71-f757-45e9-8727-7c4a9f404f94",
   "metadata": {},
   "source": [
    "The outputs should then look similar to the following: ![mgm1](notebook5_images/mgm.png)\n",
    "![mgm2](notebook5_images/mgm_2.png)\n",
    "![mgm3](notebook5_images/MGM_total.png)\n",
    "\n",
    "## Animation of the gradient-descent fitting\n",
    "\n",
    "The fitting history can also be animated.\n",
    "Set `save=True` to write the animation as a GIF."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "id": "3a79c576-3e8e-4619-b097-f863b2bf6562",
   "metadata": {},
   "outputs": [],
   "source": [
    "SA.create_animation(losses , gaussian_sums , SA.final , n_wav ,\n",
    " frame_step = 20 ,\n",
    " FPS = 15 ,\n",
    " save = False ,\n",
    " path = 'FRT00009b5a/' ,\n",
    " name = 'mgm_gradient_descent') # Animation of the fitting process"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d4b6ed6a-eeea-4210-9f3d-ae93eb11744a",
   "metadata": {},
   "source": []
  }
 ],
 "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.20"
  },
  "nbsphinx": {
   "execute": "never"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}