From cfff23f7717cb94b19eb704ac6d48f8db13d2bd2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Rapha=C3=ABl=20Robidas?= Date: Fri, 15 Dec 2023 11:11:30 -0500 Subject: [PATCH 1/2] Added an alternative way to generate first-order multiplets with better scaling --- src/nmrsim/firstorder.py | 77 ++++++++++++++++++++++++++++++++++++++-- 1 file changed, 75 insertions(+), 2 deletions(-) diff --git a/src/nmrsim/firstorder.py b/src/nmrsim/firstorder.py index 414df9a..8f00dd1 100644 --- a/src/nmrsim/firstorder.py +++ b/src/nmrsim/firstorder.py @@ -8,6 +8,7 @@ same v/J parameters that are used for second-order spin systems. See nmrsim.qm for details on these parameters. """ +from scipy.special import binom from nmrsim.math import reduce_peaks @@ -46,10 +47,10 @@ def multiplet(signal, couplings): signal : (float, float) a (frequency (Hz), intensity) tuple; couplings : [(float, int)...] - A list of (*J*, # of nuclei) tuples. The order of the tuples in + a list of (*J*, # of nuclei) tuples. The order of the tuples in couplings does not matter. e.g. to split a signal into a *dt, J* = 8, 5 Hz, use: - ``couplings = [(8, 2), (5, 3)]`` + ``couplings = [(8, 1), (5, 2)]`` Returns ------- @@ -64,6 +65,78 @@ def multiplet(signal, couplings): return sorted(reduce_peaks(res)) +def _multiplet(signal, coupling): + """ + Splits a single signal into a first-order multiplet using binomial coefficients. + + Parameters + --------- + signal : (float, float) + a (frequency (Hz), intensity) tuple; + couplings : (float, int) + a (*J*, # of nuclei) tuple. + e.g. to split a signal into a *t* with *J* = 8 Hz, use: + ``couplings = (8, 2)`` + + Returns + ------- + [(float, float)...] + a sorted peaklist for the multiplet that results from splitting the + signal by each J. + """ + + shift, intensity = signal + J, nsplit = coupling + + # Calculate the chemical shift of the first peak of the multiplet + # s -> no offset + # d -> offset by J/2 + # t -> offset by J + # ... + first_peak = shift - J * nsplit / 2 + + new_signal = [] + + # Precalculate binomial coefficients + coeffs = [binom(nsplit, i) for i in range(nsplit + 1)] + normalization_constant = intensity / sum(coeffs) + for i in range(nsplit + 1): + new_signal.append((first_peak + i * J, coeffs[i] * normalization_constant)) + + return new_signal + + +def binomial_multiplet(signal, couplings): + """ + Splits a set of signals into first-order multiplets using Pascal's triangle/binomial coefficients. + Equivalent to ```nmrsim.firstorder.multiplet``` + + Parameters + --------- + signal : (float, float) + a (frequency (Hz), intensity) tuple; + couplings : [(float, int)...] + a list of (*J*, # of nuclei) tuples. The order of the tuples in + couplings does not matter. + e.g. to split a signal into a *dt, J* = 8, 5 Hz, use: + ``couplings = [(8, 1), (5, 2)]`` + + Returns + ------- + [(float, float)...] + a sorted peaklist for the multiplet that results from splitting the + signal by each J. + """ + res = [signal] + for coupling in couplings: + new_res = [] + for sig in res: + new_res += _multiplet(sig, coupling) + res = new_res + + return sorted(reduce_peaks(res)) + + def first_order_spin_system(v, J): """ Create a first-order peaklist of several multiplets from the same v/J From 0420d94a0db44395a67982b65661355695970b90 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Rapha=C3=ABl=20Robidas?= Date: Mon, 15 Jan 2024 11:47:12 -0500 Subject: [PATCH 2/2] Replaced scipy by the standard library --- src/nmrsim/firstorder.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/nmrsim/firstorder.py b/src/nmrsim/firstorder.py index 8f00dd1..1b48f3c 100644 --- a/src/nmrsim/firstorder.py +++ b/src/nmrsim/firstorder.py @@ -8,7 +8,7 @@ same v/J parameters that are used for second-order spin systems. See nmrsim.qm for details on these parameters. """ -from scipy.special import binom +from math import comb from nmrsim.math import reduce_peaks @@ -98,7 +98,7 @@ def _multiplet(signal, coupling): new_signal = [] # Precalculate binomial coefficients - coeffs = [binom(nsplit, i) for i in range(nsplit + 1)] + coeffs = [comb(nsplit, i) for i in range(nsplit + 1)] normalization_constant = intensity / sum(coeffs) for i in range(nsplit + 1): new_signal.append((first_peak + i * J, coeffs[i] * normalization_constant))