Feat: Add support for spin-1 nuclei

.gitignore

@@ -21,4 +21,6 @@ docs/build
 docs/source/_build/
 .tox
 .DS_Store
-.coverage*
+.coverage*
+# Ignore generated NMRsim binary files
+*.npz

src/nmrsim/_classes.py

@@ -134,6 +134,10 @@ class SpinSystem:
         corresponds to the column and row order in the matrix, e.g.
         couplings[0][1] and [1]0] are the J coupling between the nuclei of
         freqs[0] and freqs[1].
+    s : array-like of float, optional (default = None)
+        A list or array containing the spin quantum number (I) for each nucleus.
+        If None, all nuclei are assumed to be spin-1/2 (I=0.5).
+        E.g., np.array([0.5, 0.5, 1.0]) for two spin-1/2 and one spin-1 nucleus.
     w : float or int (optional, default = 0.5)
         the peak width (in Hz) at half-height.
         Currently only used when SpinSystem is a component of a nmrsim.Spectrum
@@ -146,6 +150,8 @@ class SpinSystem:
     ----------
     v
     J
+    s : np.ndarray
+        Array of spin quantum numbers for each nucleus.
     w : float or int (optional, default = 0.5)
         the peak width (in Hz) at half-height.
         Currently only used when SpinSystem is a component of a nmrsim.Spectrum
@@ -160,14 +166,50 @@ class SpinSystem:
     Addition returns a Spectrum object with the SpinSystem as a component.
 
     """
+    # Add 's' to the __init__ signature with a default of None
+    def __init__(self, v, J, s=None, w=0.5, second_order=True):
 
-    def __init__(self, v, J, w=0.5, second_order=True):
         self._nuclei_number = len(v)
         self.v = v
         self.J = J
+        # Handle 's' (spin quantum numbers)
+        if s is None:
+            # Default to all spin-1/2 if 's' is not provided
+            self._s = np.full(self._nuclei_number, 0.5)
+        else:
+            # Convert to numpy array and validate
+            s_array = np.array(s)
+            if len(s_array) != self._nuclei_number:
+                raise ValueError("Length of 's' (spin quantum numbers) must match length of 'v'.")
+            if not np.all(np.isin(s_array, [0.5, 1.0])):  # Validate spin values
+                raise ValueError("Only spin quantum numbers 0.5 and 1.0 are currently supported.")
+            self._s = s_array
         self.w = w
         self._second_order = second_order
-        self._peaklist = self.peaklist()
+        self._peaklist = None
+
+    # Add a property for 's' (spin quantum numbers)
+    @property
+    def s(self):
+        """An array of the spin quantum number for each nucleus.
+
+        Returns
+        -------
+        np.ndarray
+        """
+        return self._s
+
+    # Add a setter for 's' for potential future updates, with validation
+    @s.setter
+    def s(self, s_list):
+        s_array = np.array(s_list)
+        if len(s_array) != self._nuclei_number:
+            raise ValueError("Length of 's' (spin quantum numbers) must match the number of nuclei.")
+        if not np.all(np.isin(s_array, [0.5, 1.0])):
+            raise ValueError("Only spin quantum numbers 0.5 and 1.0 are currently supported.")
+        self._s = s_array
+        # If 's' changes, the peaklist might change, so invalidate it
+        self._peaklist = None  # or self.peaklist() if you want it recomputed immediately
 
     @property
     def v(self):
@@ -186,7 +228,7 @@ def v(self, vlist):
             raise ValueError("v length must match J shape.")
         if not isinstance(vlist[0], numbers.Real):
             raise TypeError("v must be an array of numbers.")
-        self._v = vlist
+        self._v = np.array(vlist)  # Ensure it's stored as a numpy array
 
     @property
     def J(self):
@@ -213,7 +255,7 @@ def J(self, J_array):
         for i in range(m):
             if J[i, i] != 0:
                 raise ValueError("Diagonal elements of J must be 0.")
-        self._J = J_array
+        self._J = J  # Store as numpy array already
 
     @property
     def second_order(self):
@@ -236,24 +278,48 @@ def second_order(self, boolean):
         else:
             raise TypeError("second_order must be a boolean")
 
-    def peaklist(self):
+    def peaklist(self, **kwargs):
         """Return a list of (frequency, intensity) signals.
 
         Returns
         -------
         [(float, float)...]
             Array of (frequency, intensity) signals.
         """
+        if self._peaklist is not None:  # Check if peaklist is already computed
+            if kwargs:
+                pass  # Proceed to computation
+            else:
+                return self._peaklist
+
         if self._second_order:
-            return qm_spinsystem(self._v, self._J)
+            # Pass the new 's' attribute to the qm_spinsystem function
+            # Make sure qm is imported at the top of the file: from . import qm
+            # Or if it's top-level: from nmrsim import qm
+            # Pass the new 's' attribute AND **kwargs to the qm_spinsystem function
+            computed_peaklist = qm_spinsystem(self._v, self._J, s=self._s, **kwargs)
         else:
-            return first_order_spin_system(self._v, self._J)
+            # first_order_spin_system might also need 's' in the future,
+            # but for now, it's assumed to be spin-1/2 only.
+            # If first_order_spin_system also takes normalize, you'd pass kwargs here too.
+            computed_peaklist = first_order_spin_system(self._v, self._J)
+
+        # Only cache if no kwargs are passed (i.e., default unnormalized peaklist)
+        # Or, implement a more sophisticated caching strategy for normalized/cutoff variations
+        if not kwargs:
+            self._peaklist = computed_peaklist  # Cache the computed peaklist
+
+        return computed_peaklist
 
     def __eq__(self, other):
         if hasattr(other, "peaklist") and callable(other.peaklist):
+            # Ensure both peaklists are computed before comparing
             return np.allclose(self.peaklist(), other.peaklist())
+        return NotImplemented # Important for proper __eq__ behavior
 
     def __add__(self, other):
+        # Assuming Spectrum is imported correctly
+        # e.g., from nmrsim.spec import Spectrum
         if hasattr(other, "peaklist") and callable(other.peaklist):
             return Spectrum([self, other])
         else:

src/nmrsim/math.py

@@ -48,35 +48,68 @@ def add_peaks(plist):
 
 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.
+    Takes a list of (frequency, intensity) peaks and combines those whose frequencies
+    are within a specified tolerance limit.
+
+    This function is used to simplify peak lists by merging very close or
+    numerically identical peaks that should appear as a single signal in a spectrum.
 
     Parameters
-    ---------
-    plist_ : [(float, float)...]
-        A list of (x, y) tuples
+    ----------
+    plist_ : list of (float, float) tuples or np.ndarray
+        The input peak list. Can be:
+        - A standard Python list of (frequency, intensity) tuples.
+        - A 2D NumPy array where each row is a [frequency, intensity] pair.
     tolerance : float
-        tuples that differ in x by <= tolerance are combined using `add_peaks`
+        The maximum absolute difference in frequency (x-value) for two peaks
+        to be considered "close enough" to be combined. Frequencies are in Hz.
 
     Returns
     -------
-    [(float, float)...]
-        a list of (x, y) tuples where all x values differ by > `tolerance`
+    list of (float, float) tuples
+        A new list of (frequency, intensity) tuples where closely spaced
+        peaks have been combined. Frequencies are sorted in ascending order.
     """
+
+    # Convert NumPy array input to a list of tuples to ensure consistent processing
+    # by the subsequent sorting and reduction logic, which expects a list of tuples.
+    if isinstance(plist_, np.ndarray):
+        # Handle empty NumPy array explicitly to prevent errors in list comprehension.
+        if plist_.size == 0:
+            return []
+
+        # Convert each NumPy array row (e.g., [frequency, intensity]) into a tuple.
+        # This is crucial for `sorted()` to work correctly on the elements.
+        plist = [tuple(row) for row in plist_]
+    else:
+        # If the input is not a NumPy array, assume it's already a list of tuples
+        # or a compatible sequence, and proceed directly. This maintains
+        # backward compatibility with original usage.
+        plist = plist_
+
+    # Sorts the peak list by frequency (the first element of each tuple).
+    # This is essential for the reduction algorithm to correctly group
+    # adjacent peaks within the specified tolerance.
+    plist_sorted = sorted(plist)
+
     res = []
-    work = []  # an accumulator of peaks to be added
-    plist = sorted(plist_)
-    for peak in plist:
+    work = []  
+
+    for peak in plist_sorted:
         if not work:
             work.append(peak)
             continue
+
+        # Check if the current peak's frequency is within tolerance of the last
+        # peak added to the `work` group.
         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
+    # After the loop finishes, there might be peaks left in `work` that
+    # haven't been added to `res`. Combine and add them.            
+    if work:  
         res.append(add_peaks(work))
 
     return res
@@ -217,3 +250,61 @@ def get_maxima(lineshape):
             print((lineshape[0][index], val))
             res.append((lineshape[0][index], val))
     return res
+
+
+def ppm_to_hz_from_nuclei_info(ppm_positions, nuclei_types, gyromagnetic_ratios_MHz_per_T, spectrometer_1H_MHz):
+    """
+    Converts an array of ppm chemical shifts to absolute frequencies (Hz) based on
+    nuclear types and user-provided gyromagnetic ratios, given a reference 1H spectrometer frequency.
+
+    Parameters
+    ----------
+    ppm_positions : np.ndarray
+        Array of chemical shifts in ppm for each nucleus.
+    nuclei_types : list of str
+        List of nuclear symbols (e.g., '1H', '2H', '13C') corresponding to ppm_positions.
+        This list's order should match `gyromagnetic_ratios_MHz_per_T`.
+    gyromagnetic_ratios_MHz_per_T : list of float
+        List of gyromagnetic ratios in MHz/Tesla for each nucleus in the order
+        they appear in `nuclei_types`. This allows users to input exotic nuclei.
+    spectrometer_1H_MHz : float
+        The reference frequency of the spectrometer in MHz (e.g., 600.0 MHz for 1H).
+        This value, along with the hardcoded 1H gyromagnetic ratio, is used to calculate
+        the B0 magnetic field strength.
+
+    Returns
+    -------
+    np.ndarray
+        An array of absolute frequencies in Hz for each nucleus, suitable for `SpinSystem(v, ...)`.
+
+    Raises
+    ------
+    ValueError
+        If input array lengths do not match.
+    """
+
+    # Hardcoded 1H gyromagnetic ratio in MHz/Tesla
+    # This is a fundamental constant used to calculate B0 from the spectrometer's 1H frequency.
+    _GAMMA_1H_MHZ_PER_T = 42.577478461 
+
+    if not (len(ppm_positions) == len(nuclei_types) == len(gyromagnetic_ratios_MHz_per_T)):
+        raise ValueError("All input arrays (ppm_positions, nuclei_types, gyromagnetic_ratios_MHz_per_T) "
+                         "must have the same length.")
+
+    # Calculate the B0 field strength in Tesla using the hardcoded 1H gamma
+    # and the user-provided spectrometer's 1H frequency.
+    B0_Tesla = spectrometer_1H_MHz / _GAMMA_1H_MHZ_PER_T
+
+    v_calculated = []
+    for i, ppm in enumerate(ppm_positions):
+        gamma_nucleus_MHz_per_T = gyromagnetic_ratios_MHz_per_T[i]
+
+        # Calculate the Larmor frequency of *this specific nucleus* at the calculated B0 field, in Hz
+        nucleus_larmor_freq_Hz = gamma_nucleus_MHz_per_T * B0_Tesla * 1_000_000 # Convert MHz to Hz
+
+        # Calculate the actual frequency of the signal in Hz based on its ppm shift.
+        chemical_shift_in_Hz = ppm * (nucleus_larmor_freq_Hz / 1_000_000) # (Hz/ppm) * ppm
+
+        v_calculated.append(chemical_shift_in_Hz)
+
+    return np.array(v_calculated)