code/Spectra conversion script/Bad Blood 1.2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bad Blood\n",
    "\n",
    "The objective of this program is to analyse any IR spectra and extract info from their file names.\n",
    "\n",
    "Written by Mario González Jiménez  \n",
    "School of Chemistry, University of Glasgow  \n",
    "Glasgow G12 8QQ, UK.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Changelog\n",
    "\n",
    "- Version 1.0.\n",
    "    - Original Program based on ``Loco mosquito 5.0``\n",
    "- Version 1.1.\n",
    "    - Works with filenames of any number of categories.\n",
    "- Version 1.2.\n",
    "    - Solved bug processing filenames with dots"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What spectra can Bad Blood read?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Bad Blood can read any ``*.dpt`` or ``*.mzz`` file. To extract info from the name, the different categories must be separated in sections by a hyphen and **ALL** the files of the folder must have the same number of sections. If not, the program will not work. Bad Blood will interpret each section and if there are enough sections prepare a table with the number of spectra in each category."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "\n",
    "### Useful modules"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These are the modules that will be used by Bad Blood"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#we will use numpy to process the data\n",
    "import numpy as np \n",
    "# to remove the water and carbon dioxide spectra from our data\n",
    "from scipy import optimize\n",
    "\n",
    "# to make the figures\n",
    "import matplotlib.pyplot as plt \n",
    "\n",
    "# to manage the files \n",
    "import os\n",
    "# to deal with file paths on Windows, Mac and Linux\n",
    "from pathlib import Path \n",
    "# to manage the data textfiles\n",
    "import csv \n",
    "# to decompress the mzz files\n",
    "import zipfile \n",
    "import zlib\n",
    "\n",
    "# to control the program execution\n",
    "import time \n",
    "# to know the progress in the slow parts\n",
    "from tqdm import tqdm \n",
    "\n",
    "# a quick algorithm to check if all the names have the same number of sections\n",
    "def equality(listina):\n",
    "    listina = iter(listina)\n",
    "    try:\n",
    "        uno = next(listina)\n",
    "    except StopIteration:\n",
    "        return True\n",
    "    return all(uno == rest for rest in listina)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "First you must indicate the folder location. To do this, the easiest way to proceed is just going to a file in the folder where the spectra are, right-click on it, select properties, and, then, copy the text at the right of \"Location:\". **It is not necessary** to add `\\` or `/` at the end of the path."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spectra_path = input(\"Please, indicate the folder location: \")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following algorithm checks that the files in the folder are rightly named. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we include the path module that helps a lot with the / or \\ problem\n",
    "spectra_path = Path(spectra_path)\n",
    "# find all the .mzz and .dpt files in the folder (and its subfolders)\n",
    "dptfiles = []\n",
    "mzzfiles = []\n",
    "for root, dirs, files in os.walk(spectra_path):\n",
    "    for file in files:\n",
    "        if file.endswith(\".dpt\"):\n",
    "             dptfiles.append(Path(root) / Path(file))\n",
    "        elif file.endswith(\".mzz\"):\n",
    "             mzzfiles.append(Path(root) / Path(file))\n",
    "#If there are *.mzz files, we will not consider the *.dpt files\n",
    "if len(mzzfiles) > 0:\n",
    "    spectra_names = mzzfiles\n",
    "    mzzq = True\n",
    "else:\n",
    "    spectra_names = dptfiles\n",
    "    mzzq = False\n",
    "naimeision = []\n",
    "for i in spectra_names:\n",
    "    # To obtain the info from the name of the file, first we split the name in the different sections:\n",
    "    tmp = os.path.basename(i).split(\"-\")\n",
    "    tmp[-1] = tmp[-1].split(\".\")[0].split(\" \")[0]\n",
    "    tmp2 = len(tmp)\n",
    "    naimeision.append(tmp2)\n",
    "if equality(naimeision) == False:\n",
    "    nf = []\n",
    "    lf = list(set(naimeision))\n",
    "    for j in lf:\n",
    "        nf.append(naimeision.count(j))\n",
    "    jf = []\n",
    "    for k,l in enumerate(naimeision):\n",
    "        if l == lf[nf.index(min(nf))]:\n",
    "            jf.append(os.path.basename(spectra_names[k]))\n",
    "    print(\"Attention!!!!\")\n",
    "    print(\"Not all files have the same number of sections. These are probably the files that are misnamed:\")\n",
    "    print(\"\")\n",
    "    for j in jf:\n",
    "        print(\"     \"+j)\n",
    "else:\n",
    "    print(\"Everything seems all right. You may continue.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If everything went well, you can proceed with extracting the info from the names. This first part of the algorithm clasify the data and if it detects dates, let you know if there are problems with their format:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if equality(naimeision) == True:\n",
    "    tembo = [[] for jijiji in range(tmp2)]\n",
    "    for i in spectra_names:\n",
    "        # To obtain the info from the name of the file, first we split the name in the different sections:\n",
    "        tmp = os.path.basename(i).split(\"-\")\n",
    "        tmp[-1] = tmp[-1].split(\".\")[0].split(\" \")[0]\n",
    "        for k in range(len(tmp)):\n",
    "            tmp[k] = tmp[k].upper()\n",
    "        for k in range(tmp2):\n",
    "            tembo[k].append(tmp[k])\n",
    "    hakuna = [[] for jijiji in range(tmp2)]\n",
    "    matata = []\n",
    "    for i in range(len(tembo)):\n",
    "        hakuna[i] = list(set(tembo[i]))\n",
    "                \n",
    "        if len(hakuna[i]) > 0.4*len(spectra_names):\n",
    "            matata.append(\"ID\")\n",
    "        else:\n",
    "            if len(hakuna[i][0]) == 6 and (hakuna[i][0][0] == \"1\" or hakuna[i][0][0] == \"2\"):\n",
    "                matata.append(\"Dat\")\n",
    "            else:\n",
    "                matata.append(\"Cat\")\n",
    "    jf = []\n",
    "    for j, l in enumerate(matata):\n",
    "        if l == \"Dat\":\n",
    "            lf = []\n",
    "            for m,n in enumerate(tembo[j]):\n",
    "                if len(n) != 6:\n",
    "                    jf.append(os.path.basename(spectra_names[m]))\n",
    "    if len(jf) > 0:\n",
    "        print(\"Attention!!!!\")\n",
    "        print(\"There are files with wrond date format. These are probably the files that are misnamed:\")\n",
    "        print(\"\")\n",
    "        for j in jf:\n",
    "            print(\"     \"+j) \n",
    "else: \n",
    "    print(\"I told you this was not going to run unless you solve the problem with the names\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The second part of the algorithm creates the matrix with the extracted data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "matrix = []\n",
    "mbuni = [m for m,n in enumerate(matata) if n != \"Dat\"]\n",
    "kifaru = [m for m,n in enumerate(matata) if n == \"Dat\"]\n",
    "wrong_named = []\n",
    "# Now we load the spectra in a matrix\n",
    "for i in tqdm(spectra_names):\n",
    "    # To obtain the info from the name of the file, first we split the name in the different sections:\n",
    "    tmp = os.path.basename(i).split(\"-\")\n",
    "    tmp[-1] = tmp[-1].split(\".\")[0].split(\" \")[0]\n",
    "    for k in range(len(tmp)):\n",
    "        tmp[k] = tmp[k].upper()\n",
    "    # First the spectrum and its characteristics\n",
    "    if mzzq == False:\n",
    "        with open(i, 'rb') as tmp:\n",
    "            avmi = (line.replace(b'\\t',b',') for line in tmp)\n",
    "            spectrum = np.genfromtxt(avmi, delimiter=',')\n",
    "        start = spectrum[0,0]\n",
    "        end = spectrum[-1,0]\n",
    "        ls = len(spectrum)\n",
    "        spectrum = np.transpose(spectrum)[1]    \n",
    "    else:\n",
    "        with zipfile.ZipFile(i) as myzip:\n",
    "            tmpname = myzip.namelist()[0]\n",
    "            with myzip.open(tmpname) as myfile:\n",
    "                spectrum = np.genfromtxt(myfile, delimiter=',')\n",
    "        start = spectrum[0]\n",
    "        end = spectrum[1]\n",
    "        ls = int(spectrum[2])\n",
    "        spectrum = spectrum[3:]\n",
    "    # And then we incorporate all the info to the matrix\n",
    "    fisi = [[start,end,ls],spectrum] + [tmp[j] for j in mbuni]\n",
    "    if len(kifaru) == 2:\n",
    "        try:\n",
    "            colday = time.mktime(time.strptime(tmp[kifaru[0]],\"%y%m%d\"))\n",
    "        except:\n",
    "            wrong_named.append(i)\n",
    "            continue\n",
    "        try:\n",
    "            mesday = time.mktime(time.strptime(tmp[kifaru[1]],\"%y%m%d\"))\n",
    "        except:\n",
    "            wrong_named.append(i)\n",
    "            continue\n",
    "        stime = round(abs((mesday - colday) / (3600 * 24)))\n",
    "        fisi.append(stime)\n",
    "    matrix.append(fisi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is the list with the files that the algorithm detected with wrong date formats:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "wrong_named"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Number of spectra for each category"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Coming soon!!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "## Detection of spectra with low quality"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Spectra with low intensity"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If the mosquito was not well placed at the ATR's crystal, the intensity of the whole spectrum is small. Our experience says that we can use as reference the small plateau between 400 and 500 wavenumber values that the mosquito spectra usually have. Since the spectrometer with ZnSe optics only can reach to 500 wavenumber values, it doesn't matter if we extend this range to 600 cm<sup>-1</sup>. Then if the average of this reference is smaller than 0.11, the spectrum doesn't have enough quality to be scaled and, then, of course, employed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# A list of the discarted spectra will be collected:\n",
    "bad_spectra = []\n",
    "for i in range(len(matrix)):\n",
    "    # first we calculate the position of the points that comprise that section of the spectrum\n",
    "    if matrix[i][0][1] < 600 and matrix[i][0][1] > 400:\n",
    "        sta = int(round((((matrix[i][0][2] - 1) / (matrix[i][0][1] - matrix[i][0][0])) * (600 - matrix[i][0][0])) + 1)) - 1\n",
    "        end = matrix[i][0][2]\n",
    "    elif matrix[i][0][1] <= 400:\n",
    "        sta = int(round((((matrix[i][0][2] - 1) / (matrix[i][0][1] - matrix[i][0][0])) * (600 - matrix[i][0][0])) + 1)) - 1\n",
    "        end = int(round((((matrix[i][0][2] - 1) / (matrix[i][0][1] - matrix[i][0][0])) * (400 - matrix[i][0][0])) + 1)) - 1\n",
    "    else:\n",
    "        sta = 0 # if the spectrum doesn't reach 600 cm-1 we cannot prove if the spectrum has enough intensity\n",
    "        raise Exception(\"The spectrum {} doesn't reach 600 cm-1\".format(spectra_names[1]))\n",
    "    # Now we check the intensity of the spectra in that region. If is not over 0.1 we discard the spectrum\n",
    "    if np.average(matrix[i][1][sta:end]) < 0.11:\n",
    "        bad_spectra.append(\"LI: \" + str(spectra_names[i]))\n",
    "        matrix[i] = None\n",
    "if (bad_spectra) == 1:\n",
    "    print(\"1 spectrum has been discarded because its low intensity\")\n",
    "else:\n",
    "    print(str(len(bad_spectra)) + \" spectra have been discarded because their low intensity\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Spectra with abnormal background"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Bruker spectrometers sometimes measure spectra whose characteristics have been masked by a high-absorbance band covering the entire mid-infrared region. I do not know the origin of this band. I used to think it was due to the interaction of the ATR crystal with the metallic anvil, but now I'm inclined to think it's due to some vibration or shock during measurement. Sometimes these spectra pass the previous filters and it is necessary to remove them. To do it that we are going to do is to select the wavenumber with less signal from the mosquito (this wavenumber is usually 1900 cm<sup>-1</sup>) and look for outliers at that frequency."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bs = 0 # counter for the number of spectra discarderd \n",
    "# we calculate the fences of the data set based in a value we can choose (in statistics 1.5 times\n",
    "# the interquartile range is the inner fence and 3 times is the outer fence)\n",
    "l = 2.5\n",
    "# We look for the point at 1900 cm-1 and add it to the list of intensities\n",
    "li = []\n",
    "for i in range(len(matrix)):\n",
    "    if matrix[i]: #to check if we have spectra\n",
    "        # Now one would spect that the spectrum will reach 3900 so the program will not check it out.\n",
    "        sta = int(round((((matrix[i][0][2] - 1) / (matrix[i][0][1] - matrix[i][0][0])) * (1900 - matrix[i][0][0])) + 1)) - 1\n",
    "        li.append(matrix[i][1][sta])\n",
    "q3, q1 = np.percentile(li, [75 ,25])\n",
    "ir = q3 - q1\n",
    "for i in range(len(matrix)):\n",
    "    if matrix[i]: #to check if we have spectra\n",
    "        sta = int(round((((matrix[i][0][2] - 1) / (matrix[i][0][1] - matrix[i][0][0])) * (1900 - matrix[i][0][0])) + 1)) - 1\n",
    "        if matrix[i][1][sta] > (l * ir + q3) or matrix[i][1][sta] < (q1 - l * ir):\n",
    "            bs +=1\n",
    "            bad_spectra.append(\"SA: \" + str(spectra_names[i]))\n",
    "            matrix[i] = None \n",
    "if (bs) == 1:\n",
    "    print(\"1 spectrum has been discarded because it was distorted by the anvil\")\n",
    "else:\n",
    "    print(str(bs) + \" spectra have been discarded because they were distorted by the anvil\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Spectra with atmospheric interferences"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If the spectra were measured after the change of the beamsplitter or after installing the ATR and the background was not correctly measured, the spectra will be with the interference of the water and CO<sub>2</sub> spectra. In the case of water vapour, its IR spectrum has three bands with many narrow peaks that can appear as noise: one between 4000 and 3400 cm<sup>-1</sup>, other between 2200 and 1300 cm<sup>-1</sup> and the last one starts to appear below 800 cm<sup>-1</sup>. CO<sub>2</sub> has one strong band at 2345 cm<sup>-1</sup> and two smaller bands (unless you use a very fine resolution, you will not see its narrow peaks) at 3650 and 750 cm<sup>-1</sup>. Since some of this features appear in the same regions than the most interesting peaks from mosquitos, it is necessary to get rid of them. The following algorithms (originally part of a program named with another Iggy Pop song: ``Candy 0.0``) check the presence of these of these bands:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bs = 0 # counter for the number of spectra discarderd\n",
    "mycollection = []\n",
    "# Now we define a function to calculate the R-squared coefficient of the fitting of our data to a polynomial\n",
    "def rs_pf(x, y, degree):\n",
    "    coeffs = np.polyfit(x, y, degree)\n",
    "    # r-squared\n",
    "    p = np.poly1d(coeffs)\n",
    "    # fit values, and mean\n",
    "    yhat = p(x)                         # or [p(z) for z in x]\n",
    "    ybar = np.sum(y)/len(y)          # or sum(y)/len(y)\n",
    "    ssreg = np.sum((yhat-ybar)**2)   # or sum([ (yihat - ybar)**2 for yihat in yhat])\n",
    "    sstot = np.sum((y - ybar)**2)    # or sum([ (yi - ybar)**2 for yi in y])\n",
    "    results = ssreg / sstot\n",
    "\n",
    "    return results\n",
    "\n",
    "# Here take that the section of the data between 3900 and 3500 cm-1 and check if it fits well to a 5th degree polinomial\n",
    "for i in range(len(matrix)):\n",
    "    if matrix[i]: #to check if we have spectra\n",
    "        # Now one would spect that the spectrum will reach 3900 so the program will not check it out.\n",
    "        sta = int(round((((matrix[i][0][2] - 1) / (matrix[i][0][1] - matrix[i][0][0])) * (3900 - matrix[i][0][0])) + 1)) - 1\n",
    "        end = int(round((((matrix[i][0][2] - 1) / (matrix[i][0][1] - matrix[i][0][0])) * (3500 - matrix[i][0][0])) + 1)) - 1\n",
    "        # we take that data:\n",
    "        yd = matrix[i][1][sta:end]\n",
    "        xd = list(range(len(yd)))\n",
    "        rs = rs_pf(xd,yd,5)\n",
    "        # And now, if the fitting is bad, we discard the spectrum\n",
    "        if rs < 0.96:\n",
    "            bs +=1\n",
    "            bad_spectra.append(\"AI: \" + str(spectra_names[i]))\n",
    "            matrix[i] = None\n",
    "if (bs) == 1:\n",
    "    print(\"1 spectrum has been discarded because has atmospheric interferences\")\n",
    "else:\n",
    "    print(str(bs) + \" spectra have been discarded because have atmospheric interferences\") "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Number of spectra after screening"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Coming soon!!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Selection of the wavenumber values and data extraction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we proceed to finish the task extracting the intensity of the remaining spectra at the wavenumber values that we want.  We just need to indicate those wavenumber values at the next list (we have to select always more than one wavenumber). If we want to extract a range only write two numbers, the starting wavenumber and the ending wavenumber. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "wns = [3855, 3400, 3275, 2922, 2853, 1900, 1745, 1635, 1539, 1457, 1306, 1154, 1076, 1027, 880, 525]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This algorithm corrects the wavenumber values selected in case they contain wavenumber values that are not in our spectra range"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# it is very important to have the wavenumber values sorted from higher to lower\n",
    "wns.sort(reverse = True)\n",
    "\n",
    "# Now we check the lowest and highest measured wavenumber values \n",
    "a = []\n",
    "b = []\n",
    "for i in range(len(matrix)):\n",
    "    if matrix[i]:\n",
    "        a.append(matrix[i][0][0])\n",
    "        b.append(matrix[i][0][1])\n",
    "a = max(a)\n",
    "b = min(b)\n",
    "\n",
    "# If only two peaks remain, einselechta prevents the algorithm from interpreting them as a range.\n",
    "if len(wns) == 2:\n",
    "    einselechta = False\n",
    "else:\n",
    "    einselechta = True\n",
    "\n",
    "# Now we correct the wavenumber values selected that are bigger than our highest measured wavenumber\n",
    "if wns[0] > a:\n",
    "    if len(wns) == 2:\n",
    "        wns[0] = int(a)\n",
    "    else: \n",
    "        while wns[0] > a:\n",
    "            wns.pop(0)\n",
    "        if len(wns) == 1 or wns[0] < b:\n",
    "            wns.insert(0, int(a))\n",
    "    \n",
    "# And we do the same with the smaller wavenumber values:\n",
    "if wns[-1] < b:\n",
    "    if len(wns) == 2:\n",
    "        wns[-1] = int(b)\n",
    "    else:\n",
    "        while wns[-1] < b:\n",
    "            wns.pop()\n",
    "            if len(wns) == 1:\n",
    "                wns.append(int(b))   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see in a spectrum the wavenumber values selected"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We select a random spectrum:\n",
    "n = np.random.randint(0,len(matrix)-1)\n",
    "while not matrix[n]:\n",
    "    n = np.random.randint(0,len(matrix)-1)\n",
    "\n",
    "# we prepare the data\n",
    "a = matrix[n][0][0]\n",
    "b = matrix[n][0][1]\n",
    "c = matrix[n][0][2]\n",
    "xd = [a - x/c * (a-b) for x in range(c)]\n",
    "yd = matrix[n][1]\n",
    "\n",
    "# we draw the plot\n",
    "plt.figure(figsize=(14,7))\n",
    "plt.plot(xd,yd)\n",
    "plt.xlim(a, b)\n",
    "\n",
    "# and the selected wavenumber values\n",
    "if len(wns) == 2 and einselechta == False:\n",
    "    plt.axvspan(wns[0], wns[1], facecolor='papayawhip', alpha=0.5)\n",
    "else:\n",
    "    for i in wns:\n",
    "        plt.axvline(x=i, c='black', lw = 1)    \n",
    "        \n",
    "# labels\n",
    "plt.xlabel('Wavenumber')\n",
    "plt.ylabel('Absorbance')\n",
    "if len(wns) == 2 and einselechta == False:\n",
    "    plt.title('Selected Range')\n",
    "else:\n",
    "    plt.title('Selected wavenumber values')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Choose your own matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now it is the time to export the data. Since ``Bad Blood`` can deal with a lot of data, here is an algorithm to select only the spectra that we want for each category. These are the available options for each category: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "simba = [m for m,n in enumerate(matata) if n == \"Cat\"]\n",
    "for i in simba:\n",
    "    print(\"Cat \" + str(i+1) + \" - Options:\")\n",
    "    print(\"---------------- \")\n",
    "    for j in hakuna[i]:\n",
    "        print(\"      \" + j)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can make the selection by writing down in a list the options you want to use. These options have to be in lists. If you want to include all, just write ``\"all\"`` "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sel = [[\"all\"],[\"IG\", \"SU\"],[\"SF\"]]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And now, the following algorithm will extract from the spectra the desired intensities and the desired characteristics. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We start the timer\n",
    "a = time.time()\n",
    "\n",
    "# We define the variable that will contain the final data\n",
    "fida = []\n",
    "csc = 0\n",
    "ssel = 0\n",
    "\n",
    "#We convert the \"all\" statements:\n",
    "for m, n in enumerate(sel):\n",
    "    if n[0] in (\"all\",\"aLl\",\"ALL\",\"All\",\"alL\",\"aLL\"):\n",
    "        sel[m] = hakuna[simba[m]]\n",
    "\n",
    "# If we have selected a range of wavenumber values, now it is time to change it into discrete wavenumber values\n",
    "# First, we look for the boundaries of the data\n",
    "for i in range(len(matrix)):\n",
    "    if matrix[i]:\n",
    "        kk = i\n",
    "        break\n",
    "if len(wns) == 2 and einselechta == False:\n",
    "    resolution = 2\n",
    "    if matrix[kk][0][0] < wns[0]:\n",
    "        wns[0] = int(matrix[kk][0][0])\n",
    "    if matrix[kk][0][1] > wns[-1]:\n",
    "        wns[-1] = int(matrix[kk][0][-1])+1\n",
    "    wns = list(range(wns[0],wns[1],-resolution))\n",
    "\n",
    "# And start the algorithm to extract the info\n",
    "truth = True\n",
    "for i in tqdm(matrix):\n",
    "    # If that item exist\n",
    "    if i:\n",
    "        for j in range(len(sel)):\n",
    "            if truth == True and i[2+j] in sel[j]:\n",
    "                truth = True\n",
    "            else:\n",
    "                truth = False\n",
    "        if truth == True:\n",
    "            # we count the number of spectra that are removed because are too short for the range of wavenumber values selected\n",
    "            if i[0][0] >= wns[0] and i[0][1] <= wns[-1]:\n",
    "                pos = []\n",
    "                for j in wns:\n",
    "                    pos.append(int(round((((i[0][2] - 1) / (i[0][1] - i[0][0])) * (j - i[0][0])) + 1)) - 1)\n",
    "                lint = []\n",
    "                for k in pos:\n",
    "                    lint.append(i[1][k])\n",
    "                fuzz = []\n",
    "                for k in range(len(mbuni)):\n",
    "                    fuzz.append(i[2+k])\n",
    "                if len(kifaru) == 2:\n",
    "                    fuzz.append(str(int(i[2+len(mbuni)])))\n",
    "                fida.append(fuzz + lint)\n",
    "                ssel += 1\n",
    "            else:\n",
    "                csc += 1\n",
    "fida = sorted(fida)\n",
    "fluf = []\n",
    "for i in mbuni:\n",
    "    if matata[i] == \"Cat\":\n",
    "        fluf.append(matata[i] + str(i+1))\n",
    "    else:\n",
    "        fluf.append(matata[i])\n",
    "if len(kifaru) == 2:\n",
    "    fluf.append(\"StoTime\")\n",
    "fida.insert(0,fluf + wns)\n",
    "if (csc) == 0:\n",
    "    pass\n",
    "elif (csc) == 1:\n",
    "    print(\"1 spectrum has been discarded because was shorter than the selected wavenumber values\")\n",
    "else:\n",
    "    print(str(csc) + \" spectra have been discarded because were shorter than the selected wavenumber values\")  \n",
    "b = time.time()\n",
    "print(\"This last process has lasted \" + str(round(b-a,3)) + \" s. The new matrix cotains \" + str(ssel) + \" spectra.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exporting the matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we export the matrix with the info in the same folder where the data were collected"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(spectra_path.parent / Path(\"the matrix.dat\"), 'w') as file:\n",
    "    sc = csv.writer(file, delimiter='\\t')\n",
    "    for i in fida:\n",
    "            sc.writerow(i)"
   ]
  }
 ],
 "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}