Source code for nmrsim.math

"""A collection of functions for processing simulated NMR spectra.

Terms used:
signal: a pair (e.g. tuple) of frequency, intensity values
peaklist: a list (or 1D-array-like) of signals.

Provides the following functions:

* add_peaks: combines a list of signals into one signal of average frequency
  and summed intensity.

* reduce_peaks: processes a peaklist so that signals within a frequency
  tolerance are added together.

* normalize_peaklist: scales a peaklist so that intensities add to a specific
  value.
* lorentz: given a frequency, a signal and a linewidth, calculates an
  intensity. Used to calculate Lorentzian lineshapes for signals.

* get_intensity: given a lineshape and a frequency, find the intensity at the
  datapoint closest to that frequency.
"""
import numpy as np


[docs]def add_peaks(plist): """ Reduces a list of (frequency, intensity) tuples to an (average frequency, total intensity) tuple. Parameters ---------- plist: [(float, float)...] a list of (frequency, intensity) tuples Returns ------- (float, float) a tuple of (average frequency, total intensity) """ v_total = 0 i_total = 0 for v, i in plist: v_total += v i_total += i return v_total / len(plist), i_total
[docs]def reduce_peaks(plist_, tolerance=0): """ Takes a list of (x, y) tuples and adds together tuples whose values are within a certain tolerance limit. Parameters --------- plist_ : [(float, float)...] A list of (x, y) tuples tolerance : float tuples that differ in x by <= tolerance are combined using `add_peaks` Returns ------- [(float, float)...] a list of (x, y) tuples where all x values differ by > `tolerance` """ res = [] work = [] # an accumulator of peaks to be added plist = sorted(plist_) for peak in plist: if not work: work.append(peak) continue if peak[0] - work[-1][0] <= tolerance: work.append(peak) continue else: res.append(add_peaks(work)) work = [peak] if work: # process any remaining work after for loop res.append(add_peaks(work)) return res
def _normalize(intensities, n=1): """ Scale a list of intensities so that they sum to the total number of nuclei. Parameters --------- intensities : [float...] A list of intensities. n : int (optional) Number of nuclei (default = 1). """ factor = n / sum(intensities) intensities[:] = [factor * i for i in intensities]
[docs]def normalize_peaklist(peaklist, n=1): """ Normalize the intensities in a peaklist so that total intensity equals value n (nominally the number of nuclei giving rise to the signal). Parameters --------- peaklist : [(float, float)...] a list of (frequency, intensity) tuples. n : int or float (optional) total intensity to normalize to (default = 1). """ freq, int_ = [x for x, y in peaklist], [y for x, y in peaklist] _normalize(int_, n) return list(zip(freq, int_))
[docs]def lorentz(v, v0, I, w): """ A lorentz function that takes linewidth at half intensity (w) as a parameter. When `v` = `v0`, and `w` = 0.5 (Hz), the function returns intensity I. Arguments --------- v : float The frequency (x coordinate) in Hz at which to evaluate intensity (y coordinate). v0 : float The center of the distribution. I : float the relative intensity of the signal w : float the peak width at half maximum intensity Returns ------- float the intensity (y coordinate) for the Lorentzian distribution evaluated at frequency `v`. """ # Adding a height scaling factor so that peak intensities are lowered as # they are more broad. If I is the intensity with a default w of 0.5 Hz: scaling_factor = 0.5 / w # i.e. a 1 Hz wide peak will be half as high return scaling_factor * I * ((0.5 * w) ** 2 / ((0.5 * w) ** 2 + (v - v0) ** 2))
[docs]def add_lorentzians(linspace, peaklist, w): """ Given a numpy linspace, a peaklist of (frequency, intensity) tuples, and a linewidth, returns an array of y coordinates for the total line shape. Arguments --------- linspace : array-like Normally a numpy.linspace of x coordinates corresponding to frequency in Hz. peaklist : [(float, float)...] A list of (frequency, intensity) tuples. w : float Peak width at half maximum intensity. Returns ------- [float...] an array of y coordinates corresponding to intensity. """ # TODO: consider naming, and confusion with .math.add_peaks # TODO: function looks clunky. Refactor? result = lorentz(linspace, peaklist[0][0], peaklist[0][1], w) for v, i in peaklist[1:]: result += lorentz(linspace, v, i, w) return result
[docs]def get_intensity(lineshape, x): """ A crude method to find the intensity of data point closest to frequency x. Better: interpolate between two data points if match isn't exact (TODO?) Parameters ---------- lineshape : tuple of (x, y) arrays for frequency, intensity data x : frequency lookup Returns ------- float : the intensity at that frequency """ nearest_x_index = np.abs(lineshape[0] - x).argmin() return lineshape[1][nearest_x_index]
[docs]def get_maxima(lineshape): """ Crude function that returns maxima in the lineshape. Parameters ---------- lineshape : tuple of frequency, intensity arrays Returns ------- a list of (frequency, intensity) tuples for individual maxima. """ res = [] for n, val in enumerate(lineshape[1][1:-2]): index = n + 1 # start at lineshape[1][1] lastvalue = lineshape[1][index - 1] nextvalue = lineshape[1][index + 1] if lastvalue < val and nextvalue < val: print('MAXIMUM FOUND AT: ') print((lineshape[0][index], val)) res.append((lineshape[0][index], val)) return res