Improvements and bugfixes

LICENSE

@@ -1,7 +1,9 @@
 MIT License
 
-Copyright (c) 2021–2023 CHARISMA H2020 project
+Copyright (c) 2021-2024 CHARISMA H2020 project
 Copyright (c) 2022-2023 IDEAconsult Ltd.
+Copyright (c) 2024 Georgi Stefanov Georgiev
+Copyright (c) 2024 Nicolas Coca Lopez
 
 Permission is hereby granted, free of charge, to any person obtaining a copy
 of this software and associated documentation files (the "Software"), to deal

src/ramanchada2/misc/spectrum_deco/spectrum_filter.py

@@ -1,7 +1,7 @@
 #!/usr/bin/env python3
 
 from functools import wraps
-from copy import copy
+from copy import copy, deepcopy
 from ramanchada2.spectrum.spectrum import Spectrum
 from .dynamically_added import dynamically_added_filters
 import logging
@@ -13,7 +13,9 @@ def add_spectrum_filter(fun):
     @wraps(fun)
     def retf(old_spe: Spectrum, *args, **kwargs) -> Spectrum:
         new_spe = copy(old_spe)
-        new_spe._applied_processings.append(proc=fun.__name__, args=args, kwargs=kwargs)
+        new_spe._applied_processings.append(proc=fun.__name__,
+                                            args=deepcopy(args),
+                                            kwargs=deepcopy(kwargs))
         fun(old_spe, new_spe, *args, **kwargs)
         new_spe.write_cache()
         return new_spe

src/ramanchada2/misc/types/fit_peaks_result.py

@@ -8,6 +8,7 @@
 
 from ..plottable import Plottable
 from .peak_candidates import ListPeakCandidateMultiModel
+from ramanchada2.spectrum.spectrum import Spectrum
 
 
 class FitPeaksResult(UserList, Plottable):
@@ -130,3 +131,18 @@ def to_dataframe_peaks(self):
 
     def to_csv(self, path_or_buf=None, sep=',', **kwargs):
         return self.to_dataframe_peaks().to_csv(path_or_buf=path_or_buf, sep=sep, **kwargs)
+
+    def gen_fake_spectrum(self, xarr):
+        summ = np.zeros_like(xarr)
+        last_i = 0
+        last_y = 0
+        for m in self:
+            mx = m.userkws['x']
+            xi = np.searchsorted(xarr, mx)
+            evalm = m.eval()
+            summ[xi] = evalm
+            summ[np.arange(last_i, xi[0])] = np.interp(np.arange(last_i, xi[0]), [last_i, xi[0]], [last_y, evalm[0]])
+            last_i = xi[-1]
+            last_y = evalm[-1]
+        fake_spe = Spectrum(x=xarr, y=summ)
+        return fake_spe

src/ramanchada2/spectrum/baseline/add_baseline.py

@@ -14,8 +14,14 @@
 def generate_baseline(
         n_freq: int = Field(..., gt=2),
         size: int = Field(..., gt=2),
-        rng_seed: Union[int, None] = None):
-    rng = np.random.default_rng(rng_seed)
+        # validation for rng_seed is removed because
+        # it makes in-place modification impossible
+        rng_seed=None):
+    if isinstance(rng_seed, dict):
+        rng = np.random.default_rng()
+        rng.__setstate__(rng_seed)
+    else:
+        rng = np.random.default_rng(rng_seed)
     k = rng.normal(0, size, size=(2, n_freq))
     k[1][0] = 0
     z = k[0] + k[1]*1j
@@ -26,13 +32,15 @@ def generate_baseline(
     base = base[:size]
     base -= base.min()
     base /= base.max()
+    if isinstance(rng_seed, dict):
+        rng_seed.update(rng.__getstate__())
     return base
 
 
 @add_spectrum_filter
 @validate_arguments(config=dict(arbitrary_types_allowed=True))
 def add_baseline(old_spe: Spectrum, new_spe: Spectrum, /,
-                 n_freq, amplitude, pedestal, func: Union[Callable, None] = None, rng_seed=None):
+                 n_freq, amplitude, pedestal=0, func: Union[Callable, None] = None, rng_seed=None):
     """
     Add artificial baseline to the spectrum.
     A random baseline is generated in frequency domain using uniform random numbers.

src/ramanchada2/spectrum/creators/from_theoretical_lines.py

@@ -1,34 +1,40 @@
 #!/usr/bin/env python3
 
+from typing import Dict, List, Literal, Union
+
 import numpy as np
 import numpy.typing as npt
-from lmfit import Model
-from lmfit import Parameters
-from typing import Union
+from lmfit import lineshapes
 from pydantic import validate_arguments
 
-from ..spectrum import Spectrum
 from ramanchada2.misc.spectrum_deco import add_spectrum_constructor
 
+from ..spectrum import Spectrum
+
 
 @add_spectrum_constructor()
 @validate_arguments(config=dict(arbitrary_types_allowed=True))
 def from_theoretical_lines(
-        model: Model,
-        params: Parameters,
+        shapes: List[Literal[lineshapes.functions]],  # type: ignore
+        params: List[Dict],
         x: Union[int, npt.NDArray[np.float64]] = 2000):
     """
-    Generate spectrum from `lmfit` model.
+    Generate spectrum from `lmfit` shapes.
 
     Parameters
     ----------
-    model : lmfit.Model
-        the model to be used for spectrum generation
-    params : lmfit.Parameters
-        the parameters to be applied to the model
+    shapes : list of str
+        the shapes to be used for spectrum generation
+    params : list of dicts
+        shape parameters to be applied to be used with shapes
     x : Union[int, npt.NDArray[np.float64]], optional
         array with x values, by default np.array(2000)
     """
     spe = Spectrum(x=x)
-    spe.y = model.eval(params=params, x=spe.x)
+    x = spe.x
+    y = np.zeros_like(x, dtype=float)
+    for shape_name, pars in zip(shapes, params):
+        shape = getattr(lineshapes, shape_name)
+        y += shape(x=x, **pars)
+    spe.y = y
     return spe

src/ramanchada2/spectrum/filters/add_gaussian_noise.py

@@ -1,7 +1,5 @@
 #!/usr/bin/env python3
 
-from typing import Union
-
 import numpy as np
 from pydantic import validate_arguments, PositiveFloat
 
@@ -15,7 +13,9 @@ def add_gaussian_noise(
         old_spe: Spectrum,
         new_spe: Spectrum, /,
         sigma: PositiveFloat,
-        rng_seed: Union[int, None] = None):
+        # validation for rng_seed is removed because
+        # it makes in-place modification impossible
+        rng_seed=None):
     r"""
     Add gaussian noise to the spectrum.
     Random number i.i.d. $N(0, \sigma)$ is added to every sample
@@ -24,11 +24,18 @@ def add_gaussian_noise(
     ----------
     sigma : float
         sigma of the gaussian noise
-    rng_seed : int, optional
-        seed for the random generator
+    rng_seed : int or rng state, optional
+        seed for the random generator.
+        If a state is provided it is updated in-place
     """
-    rng = np.random.default_rng(rng_seed)
+    if isinstance(rng_seed, dict):
+        rng = np.random.default_rng()
+        rng.__setstate__(rng_seed)
+    else:
+        rng = np.random.default_rng(rng_seed)
     dat = old_spe.y + rng.normal(0., sigma, size=len(old_spe.y))
     if any(dat < 0):
         dat += abs(dat.min())
+    if isinstance(rng_seed, dict):
+        rng_seed.update(rng.__getstate__())
     new_spe.y = np.array(dat)

src/ramanchada2/spectrum/filters/add_gaussian_noise_drift.py

@@ -0,0 +1,74 @@
+import numpy as np
+from pydantic import PositiveFloat, confloat, validate_arguments
+
+from ramanchada2.misc.spectrum_deco import add_spectrum_filter
+
+from ..spectrum import Spectrum
+
+
+@validate_arguments(config=dict(arbitrary_types_allowed=True))
+def generate_add_gaussian_noise_drift(y, /,
+                                      sigma: PositiveFloat,
+                                      coef: confloat(ge=0, le=1),  # type: ignore [valid-type]
+                                      # validation for rng_seed is removed because
+                                      # it makes in-place modification impossible
+                                      rng_seed=None):
+    if isinstance(rng_seed, dict):
+        rng = np.random.default_rng()
+        rng.__setstate__(rng_seed)
+    else:
+        rng = np.random.default_rng(rng_seed)
+    gaus = rng.normal(0., sigma+coef/np.sqrt(2), size=len(y))
+    cs = np.cumsum(gaus)
+    # coef*sum(cs[:i]) + (1-coef)*gaus is identical to
+    # coef*sum(cs[:i-1]) + gaus
+    noise = coef*cs + gaus*(1-coef)
+    noise -= np.std(noise)
+    dat = y + noise
+    if any(dat < 0):
+        dat += abs(dat.min())
+    if isinstance(rng_seed, dict):
+        rng_seed.update(rng.__getstate__())
+    return np.array(dat)
+
+
+@add_spectrum_filter
+@validate_arguments(config=dict(arbitrary_types_allowed=True))
+def add_gaussian_noise_drift(
+        old_spe: Spectrum,
+        new_spe: Spectrum, /,
+        sigma: PositiveFloat,
+        coef: confloat(ge=0, le=1),  # type: ignore [valid-type]
+        # validation for rng_seed is removed because
+        # it makes in-place modification impossible
+        rng_seed=None):
+    r"""
+    Add cumulative gaussian noise to the spectrum.
+    Exponential-moving-average-like gaussian noise is added
+    to each sample. The goal is to mimic the low-frequency noise
+    (or random substructures in spectra).
+    The additive noise is
+    .. math::
+        a_i = coef*\sum_{j=0}^{i-1}g_j + g_i,
+    where
+    .. math::
+        g_i = \mathcal{N}(0, 1+\frac{coef}{\sqrt 2}).
+    This way drifting is possible while keeping the
+    .. math::
+        \sigma(\Delta(a)) \approx 1.
+
+    Parameters
+    ----------
+    sigma : float
+        sigma of the gaussian noise
+    coef : float in [0, 1]
+        drifting coefficient. if `coef == 0` the result is
+        identical to `add_gaussian_noise()`.
+    rng_seed : int or rng state, optional
+        seed for the random generator.
+        If a state is provided it is updated in-place
+    """
+    new_spe.y = generate_add_gaussian_noise_drift(old_spe.y,
+                                                  sigma=sigma,
+                                                  coef=coef,
+                                                  rng_seed=rng_seed)

src/ramanchada2/spectrum/filters/add_poisson_noise.py

@@ -24,10 +24,17 @@ def add_poisson_noise(
     ----------
     scale : float, optional
         scale the amplitude of the noise, by default 1
-    rng_seed : int, optional
+    rng_seed : int or rng state, optional
         seed for the random generator
+        If a state is provided it is updated in-place
     """
-    rng = np.random.default_rng(rng_seed)
+    if isinstance(rng_seed, dict):
+        rng = np.random.default_rng()
+        rng.__setstate__(rng_seed)
+    else:
+        rng = np.random.default_rng(rng_seed)
     dat = old_spe.y + [rng.normal(0., np.sqrt(i*scale)) for i in old_spe.y]
     dat[dat < 0] = 0
+    if isinstance(rng_seed, dict):
+        rng_seed.update(rng.__getstate__())
     new_spe.y = np.array(dat)

src/ramanchada2/spectrum/filters/moving_median.py

@@ -0,0 +1,42 @@
+import numpy as np
+from pydantic import PositiveInt, validate_arguments
+
+from ramanchada2.misc.spectrum_deco import add_spectrum_filter
+
+from ..spectrum import Spectrum
+
+
+@validate_arguments(config=dict(arbitrary_types_allowed=True))
+def _moving_median(s,
+                   window_size: PositiveInt = 10):
+    y = ([np.median(s[:window_size]) for i in range(window_size)] +
+         [np.median(s[i-window_size: i+window_size]) for i in range(window_size, len(s)-window_size)] +
+         [np.median(s[-window_size:]) for i in range(window_size)]
+         )
+    return np.array(y)
+
+
+@add_spectrum_filter
+@validate_arguments(config=dict(arbitrary_types_allowed=True))
+def moving_median(old_spe: Spectrum,
+                  new_spe: Spectrum, /,
+                  window_size: PositiveInt = 10):
+    """
+    Moving median filter.
+
+    Parameters
+    ----------
+    window_size : int, optional
+        by default 10
+    """
+
+    new_spe.y = _moving_median(old_spe.y, window_size)
+
+
+@add_spectrum_filter
+@validate_arguments(config=dict(arbitrary_types_allowed=True))
+def subtract_moving_median(
+        old_spe: Spectrum,
+        new_spe: Spectrum,
+        window_size: int):
+    new_spe.y = old_spe.y - _moving_median(old_spe.y, window_size)

src/ramanchada2/spectrum/peaks/find_peaks.py

@@ -80,7 +80,7 @@ def interpolate(x):
 
     peaks, props = signal.find_peaks(y_arr,
                                      prominence=prominence,
-                                     width=1,
+                                     width=width,
                                      wlen=wlen)
     peak_groups = list()
 

src/ramanchada2/spectrum/spectrum.py

@@ -10,7 +10,7 @@
 import numpy.typing as npt
 import pydantic
 from scipy.signal import convolve, savgol_coeffs, savgol_filter
-from scipy.stats import rv_histogram
+from scipy.stats import median_abs_deviation, rv_histogram
 
 from ramanchada2.io.HSDS import write_cha
 from ramanchada2.io.output.write_csv import write_csv as io_write_csv
@@ -143,6 +143,9 @@ def y(self, val: npt.NDArray[np.float64]):
     def y_noise(self):
         return self.y_noise_savgol()
 
+    def y_noise_MAD(self):
+        return median_abs_deviation(np.diff(self.y))
+
     def y_noise_savgol_DL(self, order: PositiveOddInt = PositiveOddInt(1)):
         npts = order + 2
         ydata = self.y - np.min(self.y)