add unit_conversion parameter and improve tutorial

phasorpy · cgohlke · Feb 16, 2024 · Feb 4, 2024 · Feb 5, 2024 · Feb 6, 2024
commit a6206ba9d17a23982ac388687e94478bc4a062e6
diff --git a/src/phasorpy/_phasor.pyx b/src/phasorpy/_phasor.pyx
@@ -19,7 +19,8 @@ def _phasor_from_lifetime(
     const double[::1] frequency,
     const double[:, ::1] lifetime,
     const double[:, ::1] fraction,
-    const bint is_preexp,
+    const double unit_conversion,
+    const bint preexponential,
 ):
     """Calculate phasor coordinates from lifetime components.
 
@@ -41,7 +42,10 @@ def _phasor_from_lifetime(
         Buffer of two dimensions:
             0. fractions
             1. fractions of lifetime components
-    is_preexp:
+    unit_conversion:
+        Product of `frequency` and `lifetime` units' prefix factors.
+        1e-3 for MHz and ns. 1.0 for Hz and s.
+    preexponential:
         If true, fractions are pre-exponential amplitudes, else fractional
         intensities.
 
@@ -51,7 +55,7 @@ def _phasor_from_lifetime(
         ssize_t ncomp = lifetime.shape[1]  # number lifetime components
         ssize_t ntau = lifetime.shape[0]  # number lifetimes
         ssize_t nfrac = fraction.shape[0]  # number fractions
-        double twopi = 2e-3 * M_PI
+        double twopi = 2 * M_PI * unit_conversion
         double nan = NAN  # float('NaN')
         double freq, tau, frac, sum, re, im, gs
         ssize_t f, t, s
@@ -82,7 +86,7 @@ def _phasor_from_lifetime(
                     re = 0.0
                     im = 0.0
                     sum = 0.0
-                    if is_preexp:
+                    if preexponential:
                         for s in range(ncomp):
                             sum += fraction[t, s] * lifetime[t, s]  # Fdc
                     else:
@@ -95,7 +99,7 @@ def _phasor_from_lifetime(
                     for s in range(ncomp):
                         tau = lifetime[t, s]
                         frac = fraction[t, s] / sum
-                        if is_preexp:
+                        if preexponential:
                             frac *= tau
                         tau *= freq  # omega_tau
                         gs = frac / (1.0 + tau * tau)
@@ -112,12 +116,12 @@ def _phasor_from_lifetime(
         with nogil:
             for f in range(nfreq):
                 freq = frequency[f] * twopi  # omega
-                if not is_preexp:
+                if not preexponential:
                     sum = 0.0
                     for s in range(ncomp):
                         sum += fraction[0, s]
                 for t in range(ntau):
-                    if is_preexp:
+                    if preexponential:
                         sum = 0.0
                         for s in range(ncomp):
                             sum += fraction[0, s] * lifetime[t, s]  # Fdc
@@ -130,7 +134,7 @@ def _phasor_from_lifetime(
                     for s in range(ncomp):
                         tau = lifetime[t, s]
                         frac = fraction[0, s] / sum
-                        if is_preexp:
+                        if preexponential:
                             frac *= tau
                         tau *= freq  # omega_tau
                         gs = frac / (1.0 + tau * tau)
@@ -151,7 +155,7 @@ def _phasor_from_lifetime(
                     re = 0.0
                     im = 0.0
                     sum = 0.0
-                    if is_preexp:
+                    if preexponential:
                         for s in range(ncomp):
                             sum += fraction[t, s] * lifetime[0, s]  # Fdc
                     else:
@@ -164,7 +168,7 @@ def _phasor_from_lifetime(
                     for s in range(ncomp):
                         tau = lifetime[0, s]
                         frac = fraction[t, s] / sum
-                        if is_preexp:
+                        if preexponential:
                             frac *= tau
                         tau *= freq  # omega_tau
                         gs = frac / (1.0 + tau * tau)

diff --git a/src/phasorpy/phasor.py b/src/phasorpy/phasor.py
@@ -546,7 +546,8 @@ def phasor_from_lifetime(
     /,
     fraction: ArrayLike | None = None,
     *,
-    is_preexp: bool = False,
+    preexponential: bool = False,
+    unit_conversion: float = 1e-3,
     squeeze: bool = True,
 ) -> tuple[NDArray[numpy.float64], NDArray[numpy.float64]]:
     r"""Return phasor coordinates from lifetime components.
@@ -566,9 +567,13 @@ def phasor_from_lifetime(
         Fractional intensities or pre-exponential amplitudes of the lifetime
         components. Fractions are normalized to sum to 1.
         See notes below for allowed dimensions.
-    is_preexp : bool, optional
+    preexponential : bool, optional
         If true, `fraction` values are pre-exponential amplitudes,
         else fractional intensities (default).
+    unit_conversion : float
+        Product of `frequency` and `lifetime` units' prefix factors.
+        The default is 1e-3 for MHz and ns, or Hz and ms.
+        Use 1.0 for Hz and s.
     squeeze : bool, optional
         If true (default), length-one dimensions are removed from phasor
         coordinates.
@@ -634,7 +639,8 @@ def phasor_from_lifetime(
 
     Examples
     --------
-    Phasor coordinates of a single lifetime component at 80 MHz:
+    Phasor coordinates of a single lifetime component (in ns) at a
+    frequency of 80 MHz:
 
     >>> phasor_from_lifetime(80.0, 1.9894368)  # doctest: +NUMBER
     (0.5, 0.5)
@@ -651,7 +657,7 @@ def phasor_from_lifetime(
     amplitudes:
 
     >>> phasor_from_lifetime(
-    ...     80.0, [3.9788735, 0.9947183], [0.5, 0.5], is_preexp=True
+    ...     80.0, [3.9788735, 0.9947183], [0.5, 0.5], preexponential=True
     ... )  # doctest: +NUMBER
     (0.32, 0.4)
 
@@ -678,16 +684,19 @@ def phasor_from_lifetime(
     (array([0.5, 0.5]), array([0.4, 0.5]))
 
     Phasor coordinates of multiple two-component lifetimes with specific
-    fractions at multiple frequencies:
+    fractions at multiple frequencies. Frequencies are in Hz, lifetimes in ns:
 
     >>> phasor_from_lifetime(
-    ...     [40.0, 80.0],
-    ...     [[1.0, 0.9947183], [3.9788735, 0.9947183]],
-    ...     [[0, 1], [0.5, 0.5]]
+    ...     [40e6, 80e6],
+    ...     [[1e-9, 0.9947183e-9], [3.9788735e-9, 0.9947183e-9]],
+    ...     [[0, 1], [0.5, 0.5]],
+    ...    unit_conversion=1.0
     ... )  # doctest: +NUMBER
     (array([[0.941, 0.721], [0.8, 0.5]]), array([[0.235, 0.368], [0.4, 0.4]]))
 
     """
+    if unit_conversion < 1e-16:
+        raise ValueError(f'{unit_conversion=} < 1e-16')
     frequency = numpy.array(frequency, dtype=numpy.float64, ndmin=1)
     if frequency.ndim != 1:
         raise ValueError('frequency is not one-dimensional array')
@@ -736,7 +745,9 @@ def phasor_from_lifetime(
 
     phasor = numpy.empty((2, frequency.size, nvar), dtype=numpy.float64)
 
-    _phasor_from_lifetime(phasor, frequency, lifetime, fraction, is_preexp)
+    _phasor_from_lifetime(
+        phasor, frequency, lifetime, fraction, unit_conversion, preexponential
+    )
 
     if squeeze:
         phasor = phasor.squeeze()

diff --git a/tests/test_phasor.py b/tests/test_phasor.py
@@ -493,28 +493,28 @@ def test_phasor_center_exceptions():
             ([0.72058825, 0.5], [0.36764705, 0.4]),
         ),
         # preexponential amplitudes
-        ((80.0, 0.0), {'is_preexp': True}, (1.0, 0.0)),
-        ((80.0, 1e9), {'is_preexp': True}, (0.0, 0.0)),
-        ((80.0, 3.9788735), {'is_preexp': True}, (0.2, 0.4)),
-        ((80.0, [0.0, 1e9], [0.5, 0.5]), {'is_preexp': True}, (0.0, 0.0)),
+        ((80.0, 0.0), {'preexponential': True}, (1.0, 0.0)),
+        ((80.0, 1e9), {'preexponential': True}, (0.0, 0.0)),
+        ((80.0, 3.9788735), {'preexponential': True}, (0.2, 0.4)),
+        ((80.0, [0.0, 1e9], [0.5, 0.5]), {'preexponential': True}, (0.0, 0.0)),
         (
             (80.0, [3.9788735, 0.9947183], [0.0, 1.0]),
-            {'is_preexp': True},
+            {'preexponential': True},
             (0.8, 0.4),
         ),
         (
             (80.0, [3.9788735, 0.9947183], [1.0, 0.0]),
-            {'is_preexp': True},
+            {'preexponential': True},
             (0.2, 0.4),
         ),
         (
             (80.0, [3.9788735, 0.9947183], [0.5, 0.5]),
-            {'is_preexp': True},
+            {'preexponential': True},
             (0.32, 0.4),
         ),
         (
             (80.0, [3.9788735, 0.9947183], [0.25, 0.75]),
-            {'is_preexp': True},
+            {'preexponential': True},
             (0.457143, 0.4),
         ),
         # TODO: variable lifetime, constant fraction

diff --git a/tutorials/phasorpy_phasor_from_lifetime.py b/tutorials/phasorpy_phasor_from_lifetime.py
@@ -2,8 +2,8 @@
 Phasor coordinates from lifetimes
 =================================
 
-An introduction to the :py:func:`phasorpy.phasor.phasor_from_lifetime`
-function to calculate phasor coordinates as a function of frequency,
+The :py:func:`phasorpy.phasor.phasor_from_lifetime` function is used
+to calculate phasor coordinates as a function of frequency,
 single or multiple lifetime components, and the pre-exponential amplitudes
 or fractional intensities of the components.
 
@@ -23,14 +23,14 @@
 )
 
 
-def phasor_plot(real, imag):
+def phasor_plot(real, imag, fmt='o'):
     """Plot phasor coordinates."""
     fig, ax = pyplot.subplots()
     ax.set(xlabel='real', ylabel='imag')
     ax.axis('equal')
     ax.plot(*phasor_semicircle(), color='k', lw=0.25)
     for re, im in zip(numpy.array(real, ndmin=2), numpy.array(imag, ndmin=2)):
-        ax.scatter(re, im)
+        ax.plot(re, im, fmt)
     pyplot.show()
 
 
@@ -40,10 +40,12 @@ def polar_plot(frequency, phase, modulation):
     ax.set_xscale('log', base=10)
     ax.set_xlabel('frequency (MHz)')
     ax.set_ylabel('phase (°)', color='tab:blue')
-    ax.plot(frequency, numpy.rad2deg(phase), color='tab:blue')
+    for phi in numpy.array(phase, ndmin=2).swapaxes(0, 1):
+        ax.plot(frequency, numpy.rad2deg(phi), color='tab:blue')
     ax = ax.twinx()
     ax.set_ylabel('modulation (%)', color='tab:red')
-    ax.plot(frequency, modulation * 100, color='tab:red')
+    for mod in numpy.array(modulation, ndmin=2).swapaxes(0, 1):
+        ax.plot(frequency, mod * 100, color='tab:red')
     pyplot.show()
 
 
@@ -52,9 +54,11 @@ def polar_plot(frequency, phase, modulation):
 # --------------------------
 #
 # The phasor coordinates of single-component lifetimes are located
-# on the universal circle:
+# on the universal circle.
+# For example, 3.9788735 ns and 0.9947183 ns at a frequency of 80 MHz:
 
 lifetime = numpy.array([3.9788735, 0.9947183])
+
 phasor_plot(*phasor_from_lifetime(80.0, lifetime))
 
 # %%
@@ -68,8 +72,8 @@ def polar_plot(frequency, phase, modulation):
 fraction = numpy.array(
     [[1, 0], [0.25, 0.75], [0.5, 0.5], [0.75, 0.25], [0, 1]]
 )
-phasor_plot(*phasor_from_lifetime(80.0, lifetime, fraction))
 
+phasor_plot(*phasor_from_lifetime(80.0, lifetime, fraction), fmt='o-')
 
 # %%
 # Pre-exponential amplitudes
@@ -78,45 +82,86 @@ def polar_plot(frequency, phase, modulation):
 # The phasor coordinates of two lifetime components with varying
 # pre-exponential amplitudes are also located on a line:
 
-phasor_plot(*phasor_from_lifetime(80.0, lifetime, fraction, is_preexp=True))
-
+phasor_plot(
+    *phasor_from_lifetime(80.0, lifetime, fraction, preexponential=True),
+    fmt='o-',
+)
 
 # %%
 # Lifetime distributions at multiple frequencies
 # ----------------------------------------------
 #
 # Phasor coordinates can be calculated at once for many frequencies,
-# lifetime components, and their fractions:
+# lifetime components, and their fractions.
+# As an example, lifetimes are passed in units of s and frequencies in Hz,
+# requiring to specify a unit_conversion factor:
 
-size = 100
+samples = 100
 rng = numpy.random.default_rng()
-lifetime_distribution = numpy.column_stack(
-    (
-        rng.normal(3.9788735, 0.05, size),
-        rng.normal(1.9894368, 0.05, size),
-        rng.normal(0.9947183, 0.05, size),
+lifetime_distribution = (
+    numpy.column_stack(
+        (
+            rng.normal(3.9788735, 0.05, samples),
+            rng.normal(1.9894368, 0.05, samples),
+            rng.normal(0.9947183, 0.05, samples),
+        )
     )
+    * 1e-9
 )
 fraction_distribution = numpy.column_stack(
-    (rng.random(size), rng.random(size), rng.random(size))
+    (rng.random(samples), rng.random(samples), rng.random(samples))
 )
+
 phasor_plot(
     *phasor_from_lifetime(
-        [40.0, 80.0, 160.0], lifetime_distribution, fraction_distribution
-    )
+        [40e6, 80e6, 160e6],
+        lifetime_distribution,
+        fraction_distribution,
+        unit_conversion=1.0,
+    ),
+    fmt='.',
 )
 
 # %%
-# Classic frequency-domain plots
-# ------------------------------
+# FRET efficiency
+# ---------------
 #
-# Phase-shift and demodulation of multi-component lifetimes can easily
-# be calculated as a function of the excitation light frequency:
+# The phasor coordinates of a fluorescence energy transfer donor
+# with a lifetime of 4.2 ns as a function of FRET efficiency
+# at a frequency of 80 MHz, with 90 % of the donors participating in
+# energy transfer:
+
+samples = 25
+efficiency = numpy.linspace(0.0, 1.0, samples)
+
+phasor_plot(
+    *phasor_from_lifetime(
+        80.0,
+        numpy.column_stack(
+            (
+                numpy.full(samples, 4.2),  # donor-only lifetime
+                4.2 * (1.0 - efficiency),  # donor lifetime with FRET
+            )
+        ),
+        [0.1, 0.9],
+        preexponential=True,
+    ),
+    fmt='o-',
+)
+
+# %%
+# Multi-frequency plot
+# --------------------
+#
+# Phase shift and demodulation of multi-component lifetimes can be calculated
+# as a function of the excitation light frequency and fractional intensities:
 
 frequency = numpy.logspace(-1, 4, 32)
+fraction = numpy.array([[1, 0], [0.5, 0.5], [0, 1]])
+
 polar_plot(
     frequency,
     *phasor_to_polar(
-        *phasor_from_lifetime(frequency, [3.9788735, 0.9947183], [0.8, 0.2])
+        *phasor_from_lifetime(frequency, [3.9788735, 0.9947183], fraction)
     ),
 )