Improve private parse_harmonic function (phasorpy#138)

schutyb · Jan 16, 2025 · e07c003 · e07c003
1 parent dd8ca89
commit e07c003
Show file tree

Hide file tree

Showing 7 changed files with 89 additions and 58 deletions.
diff --git a/src/phasorpy/_phasorpy.pyx b/src/phasorpy/_phasorpy.pyx
@@ -115,7 +115,7 @@ def _phasor_from_signal(
     # https://numpy.org/devdocs/reference/c-api/iterator.html
 
     if (
-        samples < 3
+        samples < 2
         or harmonics > samples // 2
         or phasor.shape[0] != harmonics * 2 + 1
         or phasor.shape[1] != signal.shape[0]

diff --git a/src/phasorpy/_utils.py b/src/phasorpy/_utils.py
@@ -249,52 +249,65 @@ def phasor_from_polar_scalar(
 
 def parse_harmonic(
     harmonic: int | Sequence[int] | Literal['all'] | str | None,
-    samples: int,
+    harmonic_max: int | None = None,
     /,
 ) -> tuple[list[int], bool]:
     """Return parsed harmonic parameter.
 
+    This function performs common, but not necessarily all, verifications
+    of user-provided `harmonic` parameter.
+
     Parameters
     ----------
     harmonic : int, list of int, 'all', or None
         Harmonic parameter to parse.
-    samples : int
-        Number of samples in signal.
-        Used to verify harmonic values and set maximum harmonic value.
+    harmonic_max : int, optional
+        Maximum value allowed in `hamonic`. Must be one or greater.
+        To verify against known number of signal samples,
+        pass ``samples // 2``.
+        If `harmonic='all'`, a range of harmonics from one to `harmonic_max`
+        (included) is returned.
 
     Returns
     -------
     harmonic : list of int
         Parsed list of harmonics.
     has_harmonic_axis : bool
-        If true, `harmonic` input parameter is an integer, else a list.
+        False if `harmonic` input parameter is a scalar integer.
 
     Raises
     ------
     IndexError
-        Any element is out of range [1..samples // 2].
+        Any element is out of range `[1..harmonic_max]`.
     ValueError
         Elements are not unique.
         Harmonic is empty.
         String input is not 'all'.
+        `harmonic_max` is smaller than 1.
     TypeError
         Any element is not an integer.
+        `harmonic` is `'all'` and `harmonic_max` is None.
 
     """
-    if samples < 3:
-        raise ValueError(f'{samples=} < 3')
+    if harmonic_max is not None and harmonic_max < 1:
+        raise ValueError(f'{harmonic_max=} < 1')
 
     if harmonic is None:
         return [1], False
 
-    harmonic_max = samples // 2
     if isinstance(harmonic, (int, numbers.Integral)):
-        if harmonic < 1 or harmonic > harmonic_max:
+        if harmonic < 1 or (
+            harmonic_max is not None and harmonic > harmonic_max
+        ):
             raise IndexError(f'{harmonic=} out of range [1..{harmonic_max}]')
         return [int(harmonic)], False
 
     if isinstance(harmonic, str):
         if harmonic == 'all':
+            if harmonic_max is None:
+                raise TypeError(
+                    f'maximum harmonic must be specified for {harmonic=!r}'
+                )
             return list(range(1, harmonic_max + 1)), True
         raise ValueError(f'{harmonic=!r} is not a valid harmonic')
 
@@ -303,10 +316,10 @@ def parse_harmonic(
         raise ValueError(f'{harmonic=} is empty')
     if h.dtype.kind not in 'iu' or h.ndim != 1:
         raise TypeError(f'{harmonic=} element not an integer')
-    if numpy.any(h < 1) or numpy.any(h > harmonic_max):
-        raise IndexError(
-            f'{harmonic=} element out of range [1..{harmonic_max}]'
-        )
+    if numpy.any(h < 1):
+        raise IndexError(f'{harmonic=} element < 1')
+    if harmonic_max is not None and numpy.any(h > harmonic_max):
+        raise IndexError(f'{harmonic=} element > {harmonic_max}]')
     if numpy.unique(h).size != h.size:
         raise ValueError(f'{harmonic=} elements must be unique')
     return h.tolist(), True

diff --git a/src/phasorpy/io.py b/src/phasorpy/io.py
@@ -290,11 +290,9 @@ def phasor_to_ometiff(
             imag = imag.reshape(1, -1)
 
     if harmonic is not None:
-        harmonic_array = numpy.atleast_1d(harmonic)
-        if harmonic_array.ndim > 1 or harmonic_array.size != nharmonic:
+        harmonic, _ = parse_harmonic(harmonic)
+        if len(harmonic) != nharmonic:
             raise ValueError('invalid harmonic')
-        samples = int(harmonic_array.max()) * 2 + 1
-        harmonic, _ = parse_harmonic(harmonic, samples)
 
     if frequency is not None:
         frequency_array = numpy.atleast_2d(frequency).astype(numpy.float64)
@@ -488,7 +486,7 @@ def phasor_from_ometiff(
 
         has_harmonic_dim = tif.series[1].ndim > tif.series[0].ndim
         nharmonics = tif.series[1].shape[0] if has_harmonic_dim else 1
-        maxharmonic = nharmonics
+        harmonic_max = nharmonics
         for i in (3, 4):
             if len(tif.series) < i + 1:
                 break
@@ -499,10 +497,10 @@ def phasor_from_ometiff(
             elif series.name == 'Phasor harmonic':
                 if not has_harmonic_dim and data.size == 1:
                     attrs['harmonic'] = int(data.item(0))
-                    maxharmonic = attrs['harmonic']
+                    harmonic_max = attrs['harmonic']
                 elif has_harmonic_dim and data.size == nharmonics:
                     attrs['harmonic'] = data.tolist()
-                    maxharmonic = max(attrs['harmonic'])
+                    harmonic_max = max(attrs['harmonic'])
                 else:
                     logger.warning(
                         f'harmonic={data} does not match phasor '
@@ -535,7 +533,7 @@ def phasor_from_ometiff(
             imag = tif.series[2].asarray()
         else:
             # specified harmonics
-            harmonic, keepdims = parse_harmonic(harmonic, 2 * maxharmonic + 1)
+            harmonic, keepdims = parse_harmonic(harmonic, harmonic_max)
             try:
                 if isinstance(harmonic_stored, list):
                     index = [harmonic_stored.index(h) for h in harmonic]
@@ -769,7 +767,7 @@ def phasor_from_simfcs_referenced(
     else:
         raise ValueError(f'file extension must be .ref or .r64, not {ext!r}')
 
-    harmonic, keep_harmonic_dim = parse_harmonic(harmonic, data.shape[0])
+    harmonic, keep_harmonic_dim = parse_harmonic(harmonic, data.shape[0] // 2)
 
     mean = data[0].copy()
     real = numpy.empty((len(harmonic),) + mean.shape, numpy.float32)

diff --git a/src/phasorpy/phasor.py b/src/phasorpy/phasor.py
@@ -104,7 +104,6 @@
 ]
 
 import math
-import numbers
 from collections.abc import Sequence
 from typing import TYPE_CHECKING
 
@@ -180,6 +179,8 @@ def phasor_from_signal(
         Else, harmonics must be at least one and no larger than half the
         number of `signal` samples along `axis`.
         The default is the first harmonic (fundamental frequency).
+        A minimum of `harmonic * 2 + 1` samples are required along `axis`
+        to calculate correct phasor coordinates at `harmonic`.
     sample_phase : array_like, optional
         Phase values (in radians) of `signal` samples along `axis`.
         If None (default), samples are assumed to be uniformly spaced along
@@ -285,7 +286,7 @@ def phasor_from_signal(
     if dtype.kind != 'f':
         raise TypeError(f'{dtype=} not supported')
 
-    harmonic, keepdims = parse_harmonic(harmonic, samples)
+    harmonic, keepdims = parse_harmonic(harmonic, samples // 2)
     num_harmonics = len(harmonic)
 
     if sample_phase is not None:
@@ -303,7 +304,7 @@ def phasor_from_signal(
         use_fft = sample_phase is None and (
             rfft is not None
             or num_harmonics > 7
-            or num_harmonics == samples // 2
+            or num_harmonics >= samples // 2
         )
 
     if use_fft:
@@ -413,7 +414,7 @@ def phasor_to_signal(
         coordinates (most commonly, lower harmonics are present if the number
         of dimensions of `mean` is one less than `real`).
         If `'all'`, the harmonics in the first axis of phasor coordinates are
-        the lower harmonics.
+        the lower harmonics necessary to synthesize `samples`.
         Else, harmonics must be at least one and no larger than half of
         `samples`.
         The phasor coordinates of missing harmonics are zeroed
@@ -462,18 +463,15 @@ def phasor_to_signal(
     array([2.2, 2.486, 0.8566, -0.4365, 0.3938])
 
     """
+    if samples < 3:
+        raise ValueError(f'{samples=} < 3')
+
     mean = numpy.array(mean, ndmin=0, copy=True)
     real = numpy.array(real, ndmin=0, copy=True)
     imag = numpy.array(imag, ndmin=1, copy=True)
 
-    if isinstance(harmonic, (int, numbers.Integral)) and harmonic == 0:
-        # harmonics are expected in the first axes of real and imag
-        samples_ = 2 * imag.shape[0]
-    else:
-        samples_ = samples
-
     harmonic_ = harmonic
-    harmonic, has_harmonic_axis = parse_harmonic(harmonic, samples_)
+    harmonic, has_harmonic_axis = parse_harmonic(harmonic, samples // 2)
 
     if real.ndim == 1 and len(harmonic) > 1 and real.shape[0] == len(harmonic):
         # single axis contains harmonic
@@ -641,7 +639,7 @@ def lifetime_to_signal(
     if harmonic is None:
         harmonic = 'all'
     all_hamonics = harmonic == 'all'
-    harmonic, _ = parse_harmonic(harmonic, samples)
+    harmonic, _ = parse_harmonic(harmonic, samples // 2)
 
     if samples < 16:
         raise ValueError(f'{samples=} < 16')

diff --git a/tests/test__utils.py b/tests/test__utils.py
@@ -130,37 +130,42 @@ def test_dilate_coordinates():
 
 def test_parse_harmonic():
     """Test parse_harmonic function."""
-    assert parse_harmonic(None, 3) == ([1], False)
-    assert parse_harmonic(1, 3) == ([1], False)
-    assert parse_harmonic(numpy.int32(1), 3) == ([1], False)
-    assert parse_harmonic([1], 3) == ([1], True)
-    assert parse_harmonic([numpy.int32(1)], 3) == (  # type: ignore[list-item]
+    assert parse_harmonic(None) == ([1], False)
+    assert parse_harmonic(None, 1) == ([1], False)
+    assert parse_harmonic(1) == ([1], False)
+    assert parse_harmonic(1, 1) == ([1], False)
+    assert parse_harmonic(numpy.int32(1), 1) == ([1], False)
+    assert parse_harmonic([1], 1) == ([1], True)
+    assert parse_harmonic([numpy.int32(1)], 1) == (  # type: ignore[list-item]
         [1],
         True,
     )
-    assert parse_harmonic([1, 2], 5) == ([1, 2], True)
-    assert parse_harmonic([2, 1], 5) == ([2, 1], True)
-    assert parse_harmonic(numpy.array([1, 2]), 5) == ([1, 2], True)
-    assert parse_harmonic('all', 5) == ([1, 2], True)
+    assert parse_harmonic([1, 2], 2) == ([1, 2], True)
+    assert parse_harmonic([2, 1], 2) == ([2, 1], True)
+    assert parse_harmonic(numpy.array([1, 2]), 2) == ([1, 2], True)
+    assert parse_harmonic('all', 1) == ([1], True)
+    assert parse_harmonic('all', 2) == ([1, 2], True)
 
     with pytest.raises(ValueError):
-        parse_harmonic(1, 2)
+        parse_harmonic(1, 0)
     with pytest.raises(IndexError):
-        parse_harmonic(0, 3)
+        parse_harmonic(0, 1)
     with pytest.raises(IndexError):
-        parse_harmonic(2, 3)
+        parse_harmonic(2, 1)
     with pytest.raises(IndexError):
-        parse_harmonic([1, 2], 3)
+        parse_harmonic([1, 2], 1)
     with pytest.raises(TypeError):
-        parse_harmonic([[1]], 3)  # type: ignore[list-item]
+        parse_harmonic([[1]], 1)  # type: ignore[list-item]
     with pytest.raises(ValueError):
-        parse_harmonic([], 3)
+        parse_harmonic([], 1)
     with pytest.raises(ValueError):
-        parse_harmonic([1, 1], 3)
+        parse_harmonic([1, 1], 1)
     with pytest.raises(ValueError):
-        parse_harmonic('alles', 3)
+        parse_harmonic('alles', 1)
     with pytest.raises(TypeError):
-        parse_harmonic(1.0, 3)
+        parse_harmonic(1.0, 1)
+    with pytest.raises(TypeError):
+        parse_harmonic('all')
 
 
 def test_chunk_iter():

diff --git a/tests/test_io.py b/tests/test_io.py
@@ -536,8 +536,8 @@ def test_phasor_to_ometiff_exceptions():
         with pytest.raises(ValueError):
             phasor_to_ometiff(filename, data[:1], data, data)
 
-        # invalid harmonic
-        with pytest.raises(ValueError):
+        # invalid harmonic, not an integer
+        with pytest.raises(TypeError):
             phasor_to_ometiff(
                 filename, *data, harmonic=[[1]]  # type: ignore[list-item]
             )

diff --git a/tests/test_phasor.py b/tests/test_phasor.py
@@ -67,11 +67,12 @@
 @pytest.mark.parametrize('use_fft', (True, False))
 def test_phasor_from_signal(use_fft):
     """Test phasor_from_signal function."""
-    sample_phase = numpy.linspace(0, 2 * math.pi, 7, endpoint=False)
+    samples = 7
+    sample_phase = numpy.linspace(0, 2 * math.pi, samples, endpoint=False)
     signal = 1.1 * (numpy.cos(sample_phase - 0.46364761) * 2 * 0.44721359 + 1)
     signal_copy = signal.copy()
 
-    # test scalar type
+    # scalar type
     mean, real, imag = phasor_from_signal(signal, use_fft=use_fft)
     assert mean.ndim == 0
     assert real.ndim == 0
@@ -204,6 +205,22 @@ def test_phasor_from_signal(use_fft):
             )
 
 
+@pytest.mark.parametrize('use_fft', (True, False))
+@pytest.mark.parametrize('samples', [2, 3])
+def test_phasor_from_signal_min_samples(samples, use_fft):
+    """Test phasor_from_signal function with two and three samples."""
+    sample_phase = numpy.linspace(0, 2 * math.pi, samples, endpoint=False)
+    signal = 1.1 * (numpy.cos(sample_phase - 0.46364761) * 2 * 0.44721359 + 1)
+
+    # a minimum of three samples is required to calculate correct 1st harmonic
+    if samples < 3:
+        with pytest.raises(ValueError):
+            phasor = phasor_from_signal(signal, use_fft=use_fft)
+    else:
+        phasor = phasor_from_signal(signal, use_fft=use_fft)
+        assert_allclose(phasor, (1.1, 0.4, 0.2), atol=1e-3)
+
+
 @pytest.mark.parametrize('use_fft', (True, False))
 @pytest.mark.parametrize(
     'shape, axis, dtype, dtype_out',