"""
phase.py
Estimate the phase of an oscillation using a waveform-based approach
"""

import numpy as np


def extrema_interpolated_phase(x, Ps, Ts, zeroxR=None, zeroxD=None):
    """
    Use peaks (phase 0) and troughs (phase pi/-pi) to estimate
    instantaneous phase. Also use rise and decay zerocrossings
    (phase -pi/2 and pi/2, respectively) if specified.

    Parameters
    ----------
    x : array-like 1d
        voltage time series
    Ps : array-like 1d
        time points of oscillatory peaks
    Ts : array-like 1d
        time points of oscillatory troughs
    zeroxR : array-like 1d
        indices at which oscillatory rising zerocrossings occur
    zeroxD : array-like 1d
        indices at which oscillatory decaying zerocrossings occur

    Returns
    -------
    pha : array-like 1d
        instantaneous phase

    Notes
    -----
    Sometimes, due to noise, extrema and zerocrossing estimation
    is poor, and for example, the same index may be assigned to
    both a peak and a decaying zerocrossing. Because of this,
    we first assign phase values by zerocrossings, and then
    may overwrite them with extrema phases.
    """

    # Initialize phase arrays
    # 2 phase arrays: trough pi and trough -pi
    L = len(x)
    t = np.arange(L)
    pha_tpi = np.zeros(L)*np.nan
    pha_tnpi = np.zeros(L)*np.nan

    # If specified, assign phases to zerocrossings
    if zeroxR is not None:
        pha_tpi[zeroxR] = -np.pi/2
        pha_tnpi[zeroxR] = -np.pi/2
    if zeroxD is not None:
        pha_tpi[zeroxD] = np.pi/2
        pha_tnpi[zeroxD] = np.pi/2

    # Define phases
    pha_tpi[Ps] = 0
    pha_tpi[Ts] = np.pi
    pha_tnpi[Ps] = 0
    pha_tnpi[Ts] = -np.pi

    # Interpolate to find all phases
    pha_tpi = np.interp(t, t[~np.isnan(pha_tpi)], pha_tpi[~np.isnan(pha_tpi)])
    pha_tnpi = np.interp(t, t[~np.isnan(pha_tnpi)], pha_tnpi[~np.isnan(pha_tnpi)])

    # For the phase time series in which the trough is negative pi:
    # Replace the decaying periods with these periods in the phase time
    # series in which the trough is pi
    diffs = np.diff(pha_tnpi)
    diffs = np.append(diffs, 99)
    pha_tnpi[diffs < 0] = pha_tpi[diffs < 0]

    return pha_tnpi

@@ -2183,16 +2183,469 @@
 s %3C 0%5D%0A%0A
+    # Assign the periods before the first empirical phase timepoint to NaN%0A    diffs = np.diff(pha_tnpi)%0A    first_empirical_idx = next(i for i, xi in enumerate(diffs) if xi %3E 0)%0A    pha_tnpi%5B:first_empirical_idx%5D = np.nan%0A%0A    # Assign the periods after the last empirical phase timepoint to NaN%0A    diffs = np.diff(pha_tnpi)%0A    last_empirical_idx = next(i for i, xi in enumerate(diffs%5B::-1%5D) if xi %3E 0)%0A    pha_tnpi%5B-last_empirical_idx+1:%5D = np.nan%0A%0A
     retu

fix code so that waveform interpolated phase has NaN edges