PeakMap

PeakMap documentation

Contributors: G. Honeybone (primary), C. Berger (PI)

Overview: the Physics of PeakMap

Before diving into any code, let’s start by understanding the relevant physics we draw upon in PeakMap. Each type of molecule has a characteristic structure – determined by which atoms it contains and how they are bonded together. This structure dictates how much energy is required to move the molecule from one energy level to another.

The field of molecular spectroscopy is dedicated to investigating these transitions, essentially by shining light through a gas and recording the wavelengths at which the light is absorbed. This effect is due to the molecule “stealing” energy from the light, and instead using it to “jump” up an energy level.

For a given molecule, the more likely a certain wavelength will jump energy levels, the more intense the transition is. In PeakMap and spectroscopy in general, we rely on visualizing this as wavenumber (inverse of wavelength) on the x-axis, and absorption intensity on the y-axis. In an unrealistic “perfect” environment, each transition would appear in this plot as a single vertical line. This is called a “stick spectrum”, and seen here:

Temperature, pressure, and other natural effects cause these lines to broaden, turning each vertical line into a combination of Gaussian and Lorentzian distributions. Spectroscopists refer to these shapes as “Voigt profiles”. Here’s an example spectrum generated by HITRAN (the database we use to create these spectra):

A molecule with a lot of physical symmetry has only a few distinct transitions, making its spectrum visually simple. See Methane doesn’t have this kind of symmetry, and its high number of rotational and vibrational degrees of freedom results in a chaotic, overlapping, and dense absorption spectrum. CO2 is an example of a “simple molecule”, where a lot of structural symmetry allows for trivially easy peak detection. Here’s a simulated CO2 spectrum:

This is where PeakMap comes in. Common peak-identification approaches such as SciPy’s find_peaks toolbox and the CLEAN algorithm often struggle to resolve peaks in methane’s dense, overlapping spectra. Instead, PeakMap uses a data-driven approach to tackle this problem. This involves a 2-stage machine learning pipeline, starting with a denoising model, and then a peak-finding model.

Now that we are familiar with the background physics, let’s dive into how PeakMap works.

PeakMap: the Code

PeakMap consists of 3 python files, all of which can be found in the github repository here in the /src folder.

File (1): data_generation.py

This file is responsible for generating the data used to train both models. We do this using the HITRAN API (HAPI), and using it to generate 5,000 spectra in the 2900-3200 cm^-1 wavenumber range, We split these spectra into 28,000 windows of 20.48 cm^-1 width including a normalized x and y axis for generalizability. We use a resolution (omega_step) of 0.01 cm^-1, and an SNR of 300. This SNR value can be increased, likely with better results. Ask Ryan Cole regarding recommended spectrometer SNR values (you could argue this SNR value to be upwards of 1000). We define a few functions in this file:

  • Generate_spectrum: this creates the individual 5000 spectra, adds noise, and normalizes the axes
  • Add_realistic_noise: includes 4 types of realistic noise to each spectra, intended to make the model more robust to real-world measurements. The noise values here can be modified, and it might be worthwhile to verify or change these.
  • Slice_and_resample: this function chops a full methane spectrum into smaller, overlapping chunks.
  • Merge_hitran_lines: merges lines that are extremely close together to their average position and intensity
  • Get_hitran_peak_positions: stores the values of each peak position according to the HITRAN line strength threshold of 1e-20. This is set by the noise floor (what we can resolve at SNR = 300), and is an important value to verify. This is one of the more arbitrary parameters in the model, and leads to decent results, but it is worth rethinking from first principles. Some food for thought: How do you define “peakiness?”. The answer to this question should motivate how we define whether a given wavenumber corresponds to a peak or not.
  • Label_chunk_from_hitran: This function is relatively unimportant, and is used mainly for diagnostics. It assigns each 20.48 cm^-1 window to one of 4 buckets based on how many peaks it contains.

The remainder of this file converts the spectra into numpy arrays, and then saves it as a compressed methane_dataset.npz file, which takes up about 3.5 gb of memory.

This file takes roughly 30 hours to run on Leavitt HPC. I would recommend running this file (and the full pipeline if curious) using a 50 spectra sample (find NUM_SAMPLES and change it from 5000 -> 50) for a 30 min run and test of the system.

For instructions on how to run this file on Leavitt, check out the Full Pipeline Reproducibility section below.

File (2): denoiser.py

This file contains the script required to train the methane spectrum denoiser model. It takes in the methane_dataset.npz file, and generates: methane_denoiser.h5 containing the weights and biases of the trained model. We define the core function: Build_unet. This compiles the model and defines its layers and shape. This is fairly optimized, but could be worth trying different architectures. This script then produces a training history plot (doesn’t natively show up when running in Leavitt), which shows the model Mean Squared Error (MSE) during each training epoch. This looks like:![][image4]

Where both training and validation loss decrease rapidly during the first few epochs, indicating successful convergence without overfitting.

File (3): peak_detector.py

[working on this currently]

Full Pipeline Reproducibility

To run the full pipeline, first make sure to install the correct dependencies. These include:

  • HAPI

  • Numpy

  • Matplotlip

  • H5py

  • Scipy

  • Tensorflow

  • Scikitlearn

  • Slurm

  • Leavitt

First step is to go to Terminal on mac and type the command:
“git clone https://github.com/ghoneybone/peakmap.git”

This clones the repo (another name for saving the folder containing the necessary files) from github and saves it to your machine. Here you have two options:

  1. Quick Pipeline Run

This takes roughly an hour or two to complete, and it is the best way to test new features or familiarize yourself with the environment. To do this, take the files from the /src folder in the PeakMap repo, and copy and paste them into jupyter hub (Ask Casey to familiarize you with metis.bates.edu jupyter hub access), vscode, or any other python development platform. Run the data generation file first, followed by the denoiser file and then the peak finder model. Be sure to wait until the previous file finishes before running the next. Change the number of samples (NUM_SAMPLES variable) in the data generation file from 5000 -> 50 or so, such that the file finishes within 45 mins. Might be worth lowering the number of training epochs for both models as well, to make the training process faster. With a quick 50 sample run, the pipeline will provide you with a strong estimate of its effectiveness, without requiring HPC access and/or a few days to converge.

  1. Full HPC Run

To fully retrain the models in the pipeline, or to get optimal results after implementing new changes, a full run of all 3 files on Leavitt HPC (or any other cluster/supercomputer) is necessary. To do this, follow these steps:

  • Open terminal, and type ssh [yourname]@leavitt.bates.edu, and type in your password when prompted
  • Now that you are inside Leavitt, type “cd /home/projects/BergerLab/[yourname]” to navigate to your personal directory. Verify with Casey that this is set up.
  • To move a file into your personal directory, type “exit” into your terminal or make sure that you are back in your local machine directory (not Leavitt). The easiest option is to locate the file you wish to move with finder, and then drag it into your terminal, giving you the full access path to the file. Then proceed to type “scp”, followed by the file path, followed by a space and “\”, and then followed by a space and “[yourname]@leavitt.bates.edu:/home/projects/BergerLab/[yourname]/
  • Hit “enter” and you have successfully copied the file into your folder in Leavitt

To run the file on leavitt, you must do so using a “slurm” script, which is a common method of securely running code on a computing cluster. My favorite way to do this is to use “vim” editor.

  • Type “vim [yourfilename].sh” and hit enter. This brings you into the editor window.
  • Type [i] to being insert mode
  • Copy and paste (or type!) your script into the window. A working script can be seen below:

#!/bin/bash
#
#SBATCH --job-name=denoiser_training # Assign a name to the job
#SBATCH --partition=defq # Which partition to use
#SBATCH --nodes=1 # Number of nodes
#SBATCH --output=testpy%j.log # Output that would print to the screen
# can be saved here
#SBATCH --error=err_testpy%j.log # Errors that arise will be saved here
pwd; hostname; date
module load anaconda/2023.03-1
conda activate myhapi
python 5_11_data_generation.py
echo “Job Finished:”
date

In the editor, it looks like this:
![][image5]

  • Once you have the script in your vim file, hit “esc” to exit insert mode, and then type “:wq” to save and quit the script.
  • Now you are out of the editor window, type in the command line “sbatch [yourfilename].sh” to submit the slurm job.
  • To view your job in the queue, type “squeue” in the command line.
  • To view file progress, type “vim testpy[yourjobnumber].log”
  • To view file errors type “vim err_testpy[yourjobnumber].log”
  • Once the file finishes running, the files will be automatically placed in your folder

Repeat this process for each of the 3 files. Historically, this will take a long time as the training for each of the files are roughly:

  • Data generation: 25 hours
  • Denoiser: 15 hours
  • Peak finder: 50 hours (may require sending the file to Casey to run, as students only have 48 hours of compute time maximum, whereas professors have ~7 days)

You have successfully run the PeakMap pipeline! Reach out to Professor Berger, or Georgie Honeybone (ghoneybone@bates.edu) with any questions. Good luck!