ronantakizawa/codereview-bench huggingface dataset:
Repository: scipy/scipy
File: scipy/signal/_spectral_py.py (Python)

Code before:
"""Tools for spectral analysis.
"""
import numpy as np
import numpy.typing as npt
from scipy import fft as sp_fft
from scipy._lib.deprecation import _deprecate_positional_args
from . import _signaltools
from ._short_time_fft import ShortTimeFFT, FFT_MODE_TYPE
from .windows import get_window
from ._arraytools import const_ext, even_ext, odd_ext, zero_ext
import warnings
from typing import cast, Literal


__all__ = ['periodogram', 'welch', 'lombscargle', 'csd', 'coherence',
           'spectrogram', 'stft', 'istft', 'check_COLA', 'check_NOLA']


@_deprecate_positional_args(version="1.19.0")
def lombscargle(
    x: npt.ArrayLike,
    y: npt.ArrayLike,
    freqs: npt.ArrayLike,
    *,
    precenter: bool | None = None,
    normalize: bool | Literal["power", "normalize", "amplitude"] = False,
    weights: npt.NDArray | None = None,
    floating_mean: bool = False,
) -> npt.NDArray:
    """
    Compute the generalized Lomb-Scargle periodogram.

Review comment (suggestion):
```suggestion
from scipy._lib.deprecation import _deprecate_positional_args, _NoValue
```

Code after:
"""Tools for spectral analysis.
"""
import numpy as np
import numpy.typing as npt
from scipy import fft as sp_fft
from scipy._lib.deprecation import _deprecate_positional_args, _NoValue
from . import _signaltools
from ._short_time_fft import ShortTimeFFT, FFT_MODE_TYPE
from .windows import get_window
from ._arraytools import const_ext, even_ext, odd_ext, zero_ext
import warnings
from typing import cast, Literal


__all__ = ['periodogram', 'welch', 'lombscargle', 'csd', 'coherence',
           'spectrogram', 'stft', 'istft', 'check_COLA', 'check_NOLA']


@_deprecate_positional_args(version="1.19.0")
def lombscargle(
    x: npt.ArrayLike,
    y: npt.ArrayLike,
    freqs: npt.ArrayLike,
    *,
    precenter: bool = _NoValue,
    normalize: bool | Literal["power", "normalize", "amplitude"] = False,
    weights: npt.NDArray | None = None,
    floating_mean: bool = False,
) -> npt.NDArray:
    """
    Compute the generalized Lomb-Scargle periodogram.
Repository: scipy/scipy
File: scipy/special/_basic.py (Python)

Code before:

    >>> import numpy as np
    >>> from scipy.special import lqn
    >>> import matplotlib.pyplot as plt

    >>> xs = np.linspace(-2, 2, 200)
    >>> n_max = 3
    >>> Qn, dQn = lqn(n_max, xs)

    Plot the Legendre functions of the second kind :math:`Q_n(x)`.

    >>> fig, ax = plt.subplots()
    >>> for n in range(n_max + 1):
    ...     ax.plot(xs, Qn[n], "-", label=rf"$n={n}$")
    >>> ax.set_xlabel(r"$x$")
    >>> ax.set_ylabel(r"$Q_n(x)$")
    >>> ax.legend()
    >>> plt.show()

    Plot the derivatives :math:`Q_n'(x)`.

    >>> fig, ax = plt.subplots()
    >>> for n in range(n_max + 1):
    ...     ax.plot(xs, dQn[n], "-", label=rf"$n={n}$")
    >>> ax.set_xlabel(r"$x$")
    >>> ax.set_ylabel(r"$Q_n'(x)$")
    >>> ax.legend()
    >>> plt.show()
    """
    n = _nonneg_int_or_fail(n, 'n', strict=False)
    if (n < 1):
        n1 = 1
    else:
        n1 = n

    z = np.asarray(z)
    if (not np.issubdtype(z.dtype, np.inexact)):
        z = z.astype(float)

    if np.iscomplexobj(z):
        qn = np.empty((n1 + 1,) + z.shape, dtype=np.complex128)
    else:
        qn = np.empty((n1 + 1,) + z.shape, dtype=np.float64)
    qd = np.empty_like(qn)
    if (z.ndim == 0):
        _lqn(z, out=(qn, qd))
    else:
          # new axes must be last for the ufunc
        _lqn(z,
             out=(np.moveaxis(qn, 0, -1),
                  np.moveaxis(qd, 0, -1)))

Review comment (suggestion):
This can be made a little more compact
```suggestion
    >>> fig, ax = plt.subplots()
    >>> ax.plot(xs, dQn[n], "-")
    >>> ax.set_xlabel(r"$x$")
    >>> ax.set_ylabel(r"$Q_n'(x)$")
    >>> ax.legend([fr"$n={n}$" for n in range(n_max + 1)])
    >>> plt.show()
```

Code after:

    >>> import numpy as np
    >>> from scipy.special import lqn
    >>> import matplotlib.pyplot as plt

    >>> xs = np.linspace(-2, 2, 200)
    >>> n_max = 3
    >>> Qn, dQn = lqn(n_max, xs)

    Plot the Legendre functions of the second kind :math:`Q_n(x)`.

    >>> fig, ax = plt.subplots()
    >>> ax.plot(xs, Qn.T, "-")
    >>> ax.set_xlabel(r"$x$")
    >>> ax.set_ylabel(r"$Q_n(x)$")
    >>> ax.legend([fr"$n={n}$" for n in range(n_max + 1)])
    >>> plt.show()

    Plot the derivatives :math:`Q_n'(x)`.

    >>> fig, ax = plt.subplots()
    >>> ax.plot(xs, dQn.T, "-")
    >>> ax.set_xlabel(r"$x$")
    >>> ax.set_ylabel(r"$Q_n'(x)$")
    >>> ax.legend([fr"$n={n}$" for n in range(n_max + 1)])
    >>> plt.show()
    """
    n = _nonneg_int_or_fail(n, 'n', strict=False)
    if (n < 1):
        n1 = 1
    else:
        n1 = n

    z = np.asarray(z)
    if (not np.issubdtype(z.dtype, np.inexact)):
        z = z.astype(float)

    if np.iscomplexobj(z):
        qn = np.empty((n1 + 1,) + z.shape, dtype=np.complex128)
    else:
        qn = np.empty((n1 + 1,) + z.shape, dtype=np.float64)
    qd = np.empty_like(qn)
    if (z.ndim == 0):
        _lqn(z, out=(qn, qd))
    else:
          # new axes must be last for the ufunc
        _lqn(z,
             out=(np.moveaxis(qn, 0, -1),
                  np.moveaxis(qd, 0, -1)))

    return qn[:(n+1)], qd[:(n+1)]
Repository: scipy/scipy
File: scipy/special/_basic.py (Python)

Code before:

    Examples
    --------
    Compute :math:`Q_n(x)` and its derivatives on an interval.

    >>> import numpy as np
    >>> from scipy.special import lqn
    >>> import matplotlib.pyplot as plt

    >>> xs = np.linspace(-2, 2, 200)
    >>> n_max = 3
    >>> Qn, dQn = lqn(n_max, xs)

    Plot the Legendre functions of the second kind :math:`Q_n(x)`.

    >>> fig, ax = plt.subplots()
    >>> ax.plot(xs, Qn.T, "-", label=rf"$n={n}$")
    >>> ax.set_xlabel(r"$x$")
    >>> ax.set_ylabel(r"$Q_n(x)$")
    >>> ax.legend([fr"$n={n}$" for n in range(n_max + 1)])
    >>> plt.show()

    Plot the derivatives :math:`Q_n'(x)`.

    >>> fig, ax = plt.subplots()
    >>> ax.plot(xs, dQn.T, "-", label=rf"$n={n}$")
    >>> ax.set_xlabel(r"$x$")
    >>> ax.set_ylabel(r"$Q_n'(x)$")
    >>> ax.legend([fr"$n={n}$" for n in range(n_max + 1)])
    >>> plt.show()
    """
    n = _nonneg_int_or_fail(n, 'n', strict=False)
    if (n < 1):
        n1 = 1
    else:
        n1 = n

    z = np.asarray(z)
    if (not np.issubdtype(z.dtype, np.inexact)):
        z = z.astype(float)

    if np.iscomplexobj(z):
        qn = np.empty((n1 + 1,) + z.shape, dtype=np.complex128)
    else:
        qn = np.empty((n1 + 1,) + z.shape, dtype=np.float64)
    qd = np.empty_like(qn)
    if (z.ndim == 0):
        _lqn(z, out=(qn, qd))
    else:
          # new axes must be last for the ufunc
        _lqn(z,

Review comment (suggestion):
```suggestion
    >>> ax.plot(xs, Qn.T, "-")
    >>> ax.set_xlabel(r"$x$")
    >>> ax.set_ylabel(r"$Q_n(x)$")
    >>> ax.legend([fr"$n={n}$" for n in range(n_max + 1)])
    >>> plt.show()

    Plot the derivatives :math:`Q_n'(x)`.

    >>> fig, ax = plt.subplots()
    >>> ax.plot(xs, dQn.T, "-")
```

Code after:

    Examples
    --------
    Compute :math:`Q_n(x)` and its derivatives on an interval.

    >>> import numpy as np
    >>> from scipy.special import lqn
    >>> import matplotlib.pyplot as plt

    >>> xs = np.linspace(-2, 2, 200)
    >>> n_max = 3
    >>> Qn, dQn = lqn(n_max, xs)

    Plot the Legendre functions of the second kind :math:`Q_n(x)`.

    >>> fig, ax = plt.subplots()
    >>> ax.plot(xs, Qn.T, "-")
    >>> ax.set_xlabel(r"$x$")
    >>> ax.set_ylabel(r"$Q_n(x)$")
    >>> ax.legend([fr"$n={n}$" for n in range(n_max + 1)])
    >>> plt.show()

    Plot the derivatives :math:`Q_n'(x)`.

    >>> fig, ax = plt.subplots()
    >>> ax.plot(xs, dQn.T, "-")
    >>> ax.set_xlabel(r"$x$")
    >>> ax.set_ylabel(r"$Q_n'(x)$")
    >>> ax.legend([fr"$n={n}$" for n in range(n_max + 1)])
    >>> plt.show()
    """
    n = _nonneg_int_or_fail(n, 'n', strict=False)
    if (n < 1):
        n1 = 1
    else:
        n1 = n

    z = np.asarray(z)
    if (not np.issubdtype(z.dtype, np.inexact)):
        z = z.astype(float)

    if np.iscomplexobj(z):
        qn = np.empty((n1 + 1,) + z.shape, dtype=np.complex128)
    else:
        qn = np.empty((n1 + 1,) + z.shape, dtype=np.float64)
    qd = np.empty_like(qn)
    if (z.ndim == 0):
        _lqn(z, out=(qn, qd))
    else:
          # new axes must be last for the ufunc
        _lqn(z,
Repository: scipy/scipy
File: scipy/special/_add_newdocs.py (Python)

Code before:
    out : ndarray, optional
        Optional output array for the function results.

    Returns
    -------
    pmv : scalar or ndarray
        Value of the associated Legendre function.

    Notes
    -----
    Note that this implementation includes the Condon-Shortley phase.

    References
    ----------
    .. [1] S. Zhang and J.M. Jin, "Computation of Special Functions", Wiley 1996.

    Examples
    --------
    Compute a single associated Legendre function value:

    >>> from scipy.special import lpmv
    >>> lpmv(1, 2.5, 0.3)
    -0.1292299

    Plot the associated Legendre functions for orders ``m = 1`` and
    ``m = 2`` and degrees ``v = 0, \dots, 5``:

    >>> import numpy as np
    >>> import matplotlib.pyplot as plt
    >>> nus = np.arange(0, 6)
    >>> xs = np.linspace(-1+1e-5, 1-1e-5, 101)

    Order 1.

    >>> fig, ax = plt.subplots()
    >>> ax.plot(xs, np.array([lpmv(1, nu, xs) for nu in nus]).T, "-")
    >>> ax.legend([rf"$\nu={nu}$" for nu in nus])
    >>> ax.set_xlabel(r"$x$")
    >>> ax.set_ylabel(r"$P_\nu^1(x)$")
    >>> plt.show()

    Order 2.

    >>> fig, ax = plt.subplots()
    >>> ax.plot(xs, np.array([lpmv(2, nu, xs) for nu in nus]).T, "-")
    >>> ax.legend([rf"$\nu={nu}$" for nu in nus])
    >>> ax.set_xlabel(r"$x$")
    >>> ax.set_ylabel(r"$P_\nu^2(x)$")
    >>> plt.show()
    """)


Review comment (suggestion):
```suggestion
    ``m = 2`` and degrees ``v = 0, ..., 5``:
```

Code after:
    out : ndarray, optional
        Optional output array for the function results.

    Returns
    -------
    pmv : scalar or ndarray
        Value of the associated Legendre function.

    Notes
    -----
    Note that this implementation includes the Condon-Shortley phase.

    References
    ----------
    .. [1] S. Zhang and J.M. Jin, "Computation of Special Functions", Wiley 1996.

    Examples
    --------
    Compute a single associated Legendre function value:

    >>> from scipy.special import lpmv
    >>> lpmv(1, 2.5, 0.3)
    -0.1292299

    Plot the associated Legendre functions for orders ``m = 1`` and
    ``m = 2`` and degrees ``v = 0, ..., 5``:

    >>> import numpy as np
    >>> import matplotlib.pyplot as plt
    >>> nus = np.arange(0, 6)
    >>> xs = np.linspace(-1+1e-5, 1-1e-5, 101)

    Order 1.

    >>> fig, ax = plt.subplots()
    >>> ax.plot(xs, lpmv(1, nus, xs[:, np.newaxis]), "-")
    >>> ax.legend([rf"$\nu={nu}$" for nu in nus])
    >>> ax.set_xlabel(r"$x$")
    >>> ax.set_ylabel(r"$P_\nu^1(x)$")
    >>> plt.show()

    Order 2.

    >>> fig, ax = plt.subplots()
    >>> ax.plot(xs, lpmv(2, nus, xs[:, np.newaxis]), "-")
    >>> ax.legend([rf"$\nu={nu}$" for nu in nus])
    >>> ax.set_xlabel(r"$x$")
    >>> ax.set_ylabel(r"$P_\nu^2(x)$")
    >>> plt.show()
    """)
Repository: scipy/scipy
File: scipy/special/_add_newdocs.py (Python)

Code before:
    Note that this implementation includes the Condon-Shortley phase.

    References
    ----------
    .. [1] S. Zhang and J.M. Jin, "Computation of Special Functions", Wiley 1996.

    Examples
    --------
    Compute a single associated Legendre function value:

    >>> from scipy.special import lpmv
    >>> lpmv(1, 2.5, 0.3)
    -0.1292299

    Plot the associated Legendre functions for orders ``m = 1`` and
    ``m = 2`` and degrees ``v = 0, \dots, 5``:

    >>> import numpy as np
    >>> import matplotlib.pyplot as plt
    >>> nus = np.arange(0, 6)
    >>> xs = np.linspace(-1+1e-5, 1-1e-5, 101)

    Order 1.

    >>> fig, ax = plt.subplots()
    >>> ax.plot(xs, np.array([lpmv(1, nu, xs) for nu in nus]).T, "-")
    >>> ax.legend([rf"$\nu={nu}$" for nu in nus])
    >>> ax.set_xlabel(r"$x$")
    >>> ax.set_ylabel(r"$P_\nu^1(x)$")
    >>> plt.show()

    Order 2.

    >>> fig, ax = plt.subplots()
    >>> ax.plot(xs, np.array([lpmv(2, nu, xs) for nu in nus]).T, "-")
    >>> ax.legend([rf"$\nu={nu}$" for nu in nus])
    >>> ax.set_xlabel(r"$x$")
    >>> ax.set_ylabel(r"$P_\nu^2(x)$")
    >>> plt.show()
    """)

add_newdoc("nbdtr",
    r"""
    nbdtr(k, n, p, out=None)

    Negative binomial cumulative distribution function.

    Returns the sum of the terms 0 through `k` of the negative binomial
    distribution probability mass function,

    .. math::

Review comment (suggestion):
We can use broadcasting rather than a list comprehension
```suggestion
    >>> ax.plot(xs, lpmv(1, nu, xs[:, np.newaxis]), "-")
```

Code after:
    Note that this implementation includes the Condon-Shortley phase.

    References
    ----------
    .. [1] S. Zhang and J.M. Jin, "Computation of Special Functions", Wiley 1996.

    Examples
    --------
    Compute a single associated Legendre function value:

    >>> from scipy.special import lpmv
    >>> lpmv(1, 2.5, 0.3)
    -0.1292299

    Plot the associated Legendre functions for orders ``m = 1`` and
    ``m = 2`` and degrees ``v = 0, ..., 5``:

    >>> import numpy as np
    >>> import matplotlib.pyplot as plt
    >>> nus = np.arange(0, 6)
    >>> xs = np.linspace(-1+1e-5, 1-1e-5, 101)

    Order 1.

    >>> fig, ax = plt.subplots()
    >>> ax.plot(xs, lpmv(1, nus, xs[:, np.newaxis]), "-")
    >>> ax.legend([rf"$\nu={nu}$" for nu in nus])
    >>> ax.set_xlabel(r"$x$")
    >>> ax.set_ylabel(r"$P_\nu^1(x)$")
    >>> plt.show()

    Order 2.

    >>> fig, ax = plt.subplots()
    >>> ax.plot(xs, lpmv(2, nus, xs[:, np.newaxis]), "-")
    >>> ax.legend([rf"$\nu={nu}$" for nu in nus])
    >>> ax.set_xlabel(r"$x$")
    >>> ax.set_ylabel(r"$P_\nu^2(x)$")
    >>> plt.show()
    """)

add_newdoc("nbdtr",
    r"""
    nbdtr(k, n, p, out=None)

    Negative binomial cumulative distribution function.

    Returns the sum of the terms 0 through `k` of the negative binomial
    distribution probability mass function,

    .. math::
Repository: scipy/scipy
File: scipy/special/_add_newdocs.py (Python)

Code before:

    >>> from scipy.special import lpmv
    >>> lpmv(1, 2.5, 0.3)
    -0.1292299

    Plot the associated Legendre functions for orders ``m = 1`` and
    ``m = 2`` and degrees ``v = 0, \dots, 5``:

    >>> import numpy as np
    >>> import matplotlib.pyplot as plt
    >>> nus = np.arange(0, 6)
    >>> xs = np.linspace(-1+1e-5, 1-1e-5, 101)

    Order 1.

    >>> fig, ax = plt.subplots()
    >>> ax.plot(xs, np.array([lpmv(1, nu, xs) for nu in nus]).T, "-")
    >>> ax.legend([rf"$\nu={nu}$" for nu in nus])
    >>> ax.set_xlabel(r"$x$")
    >>> ax.set_ylabel(r"$P_\nu^1(x)$")
    >>> plt.show()

    Order 2.

    >>> fig, ax = plt.subplots()
    >>> ax.plot(xs, np.array([lpmv(2, nu, xs) for nu in nus]).T, "-")
    >>> ax.legend([rf"$\nu={nu}$" for nu in nus])
    >>> ax.set_xlabel(r"$x$")
    >>> ax.set_ylabel(r"$P_\nu^2(x)$")
    >>> plt.show()
    """)

add_newdoc("nbdtr",
    r"""
    nbdtr(k, n, p, out=None)

    Negative binomial cumulative distribution function.

    Returns the sum of the terms 0 through `k` of the negative binomial
    distribution probability mass function,

    .. math::

        F = \sum_{j=0}^k {{n + j - 1}\choose{j}} p^n (1 - p)^j.

    In a sequence of Bernoulli trials with individual success probabilities
    `p`, this is the probability that `k` or fewer failures precede the nth
    success.

    Parameters
    ----------

Review comment (suggestion):
```suggestion
    >>> ax.plot(xs,lpmv(2, nu, xs[:, np.newaxis]), "-")
```

Code after:

    >>> from scipy.special import lpmv
    >>> lpmv(1, 2.5, 0.3)
    -0.1292299

    Plot the associated Legendre functions for orders ``m = 1`` and
    ``m = 2`` and degrees ``v = 0, ..., 5``:

    >>> import numpy as np
    >>> import matplotlib.pyplot as plt
    >>> nus = np.arange(0, 6)
    >>> xs = np.linspace(-1+1e-5, 1-1e-5, 101)

    Order 1.

    >>> fig, ax = plt.subplots()
    >>> ax.plot(xs, lpmv(1, nus, xs[:, np.newaxis]), "-")
    >>> ax.legend([rf"$\nu={nu}$" for nu in nus])
    >>> ax.set_xlabel(r"$x$")
    >>> ax.set_ylabel(r"$P_\nu^1(x)$")
    >>> plt.show()

    Order 2.

    >>> fig, ax = plt.subplots()
    >>> ax.plot(xs, lpmv(2, nus, xs[:, np.newaxis]), "-")
    >>> ax.legend([rf"$\nu={nu}$" for nu in nus])
    >>> ax.set_xlabel(r"$x$")
    >>> ax.set_ylabel(r"$P_\nu^2(x)$")
    >>> plt.show()
    """)

add_newdoc("nbdtr",
    r"""
    nbdtr(k, n, p, out=None)

    Negative binomial cumulative distribution function.

    Returns the sum of the terms 0 through `k` of the negative binomial
    distribution probability mass function,

    .. math::

        F = \sum_{j=0}^k {{n + j - 1}\choose{j}} p^n (1 - p)^j.

    In a sequence of Bernoulli trials with individual success probabilities
    `p`, this is the probability that `k` or fewer failures precede the nth
    success.

    Parameters
    ----------
Repository: scipy/scipy
File: scipy/_lib/_util.py (Python)

Code before:
            results = list(zip(*results))

            # Reshape results
            for i, result in enumerate(results):
                result = xp.stack(result)
                core_shape = result.shape[1:]
                results[i] = xp.reshape(result, batch_shape + core_shape)

            # Assume `result` should be a single array if there is only one element or
            # a `tuple` otherwise. This is easily generalized by allowing the
            # contributor to pass an `pack_result` callable to the decorator factory.
            return results[0] if len(results) == 1 else results

        doc = FunctionDoc(wrapper)
        doc['Extended Summary'].append(_batch_note.rstrip())
        wrapper.__doc__ = str(doc).split("\n", 1)[1].lstrip(" \n")  # remove signature

        return wrapper
    return decorator


def np_vecdot(x1, x2, /, *, axis=-1):
    # `np.vecdot` has advantages (e.g. see gh-22462), so let's use it when
    # available. As functions are translated to Array API, `np_vecdot` can be
    # replaced with `xp.vecdot`.
    return np.vecdot(x1, x2, axis=axis)


def _dedent_for_py313(s):
    """Apply textwrap.dedent to s for Python versions 3.13 or later."""
    return s if sys.version_info < (3, 13) else textwrap.dedent(s)


def broadcastable(shape_a: tuple[int, ...], shape_b: tuple[int, ...]) -> bool:
    """Check if two shapes are broadcastable."""
    return all(
        (m == n) or (m == 1) or (n == 1) for m, n in zip(shape_a[::-1], shape_b[::-1])
    )


Review comment (suggestion):
I think we should probably just delete `np_vecdot` and change all uses to `np.vecdot`

Code after:
            results = list(zip(*results))

            # Reshape results
            for i, result in enumerate(results):
                result = xp.stack(result)
                core_shape = result.shape[1:]
                results[i] = xp.reshape(result, batch_shape + core_shape)

            # Assume `result` should be a single array if there is only one element or
            # a `tuple` otherwise. This is easily generalized by allowing the
            # contributor to pass an `pack_result` callable to the decorator factory.
            return results[0] if len(results) == 1 else results

        doc = FunctionDoc(wrapper)
        doc['Extended Summary'].append(_batch_note.rstrip())
        wrapper.__doc__ = str(doc).split("\n", 1)[1].lstrip(" \n")  # remove signature

        return wrapper
    return decorator


def _dedent_for_py313(s):
    """Apply textwrap.dedent to s for Python versions 3.13 or later."""
    return s if sys.version_info < (3, 13) else textwrap.dedent(s)


def broadcastable(shape_a: tuple[int, ...], shape_b: tuple[int, ...]) -> bool:
    """Check if two shapes are broadcastable."""
    return all(
        (m == n) or (m == 1) or (n == 1) for m, n in zip(shape_a[::-1], shape_b[::-1])
    )
Repository: scipy/scipy
File: scipy/signal/_waveforms.py (Python)

Code before:
    >>> plt.ylim(-2, 2)

    A pulse-width modulated sine wave:

    >>> plt.figure()
    >>> sig = np.sin(2 * np.pi * t)
    >>> pwm = signal.square(2 * np.pi * 30 * t, duty=(sig + 1)/2)
    >>> plt.subplot(2, 1, 1)
    >>> plt.plot(t, sig)
    >>> plt.subplot(2, 1, 2)
    >>> plt.plot(t, pwm)
    >>> plt.ylim(-1.5, 1.5)

    """
    xp = array_namespace(t, duty)
    t, w = xp_promote(t, duty, xp=xp)

    y = xp.zeros(t.shape, dtype=xp.float64)

    # width must be between 0 and 1 inclusive
    mask1 = (w > 1) | (w < 0)
    y = xpx.at(y, mask1).set(xp.nan)

    # on the interval 0 to duty*2*pi function is 1
    tmod = t % (2 * xp.pi)
    mask2 = ~mask1 & (tmod < w * 2 * pi)
    y = xpx.at(y, mask2).set(1)

    # on the interval duty*2*pi to 2*pi function is
    #  (pi*(w+1)-tmod) / (pi*(1-w))
    mask3 = ~mask1 & ~mask2
    y = xpx.at(y, mask3).set(-1)
    return y


def gausspulse(t, fc=1000, bw=0.5, bwr=-6, tpr=-60, retquad=False,
               retenv=False):
    """
    Return a Gaussian modulated sinusoid:

        ``exp(-a t^2) exp(1j*2*pi*fc*t).``

    If `retquad` is True, then return the real and imaginary parts
    (in-phase and quadrature).
    If `retenv` is True, then return the envelope (unmodulated signal).
    Otherwise, return the real part of the modulated sinusoid.

    Parameters
    ----------
    t : ndarray or the string 'cutoff'
        Input array.

Review comment (suggestion):
```suggestion
    mask2 = ~mask1 & (tmod < w*2*xp.pi)
```

Code after:
    >>> plt.ylim(-2, 2)

    A pulse-width modulated sine wave:

    >>> plt.figure()
    >>> sig = np.sin(2 * np.pi * t)
    >>> pwm = signal.square(2 * np.pi * 30 * t, duty=(sig + 1)/2)
    >>> plt.subplot(2, 1, 1)
    >>> plt.plot(t, sig)
    >>> plt.subplot(2, 1, 2)
    >>> plt.plot(t, pwm)
    >>> plt.ylim(-1.5, 1.5)

    """
    xp = array_namespace(t, duty)
    t, w = xp_promote(t, duty, xp=xp, force_floating=True, broadcast=True)

    y = xp.zeros(t.shape, dtype=t.dtype)

    # width must be between 0 and 1 inclusive
    mask1 = (w > 1) | (w < 0)
    y = xpx.at(y, mask1).set(xp.nan)

    # on the interval 0 to duty*2*pi function is 1
    tmod = t % (2 * xp.pi)
    mask2 = ~mask1 & (tmod < w*2*xp.pi)
    y = xpx.at(y, mask2).set(1)

    # on the interval duty*2*pi to 2*pi function is -1
    mask3 = ~(mask1 | mask2)
    y = xpx.at(y, mask3).set(-1)
    return y


def gausspulse(t, fc=1000, bw=0.5, bwr=-6, tpr=-60, retquad=False,
               retenv=False):
    """
    Return a Gaussian modulated sinusoid:

        ``exp(-a t^2) exp(1j*2*pi*fc*t).``

    If `retquad` is True, then return the real and imaginary parts
    (in-phase and quadrature).
    If `retenv` is True, then return the envelope (unmodulated signal).
    Otherwise, return the real part of the modulated sinusoid.

    Parameters
    ----------
    t : ndarray or the string 'cutoff'
        Input array.
    fc : float, optional