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

Code before:
"""
Solve the orthogonal Procrustes problem.

"""

from scipy._lib._util import _apply_over_batch
from ._decomp_svd import svd
from scipy._lib._array_api import array_namespace, xp_capabilities, _asarray, is_numpy


__all__ = ['orthogonal_procrustes']


@xp_capabilities(
    jax_jit=False,
    skip_backends=[("dask.array", "full_matrics=True is not supported by dask")],
)
@_apply_over_batch(('A', 2), ('B', 2))
def orthogonal_procrustes(A, B, check_finite=True):
    """
    Compute the matrix solution of the orthogonal (or unitary) Procrustes problem.

    Given matrices `A` and `B` of the same shape, find an orthogonal (or unitary in
    the case of complex input) matrix `R` that most closely maps `A` to `B` using the
    algorithm given in [1]_.

    Parameters
    ----------
    A : (M, N) array_like
        Matrix to be mapped.
    B : (M, N) array_like
        Target matrix.
    check_finite : bool, optional
        Whether to check that the input matrices contain only finite numbers.
        Disabling may give a performance gain, but may result in problems
        (crashes, non-termination) if the inputs do contain infinities or NaNs.

    Returns
    -------
    R : (N, N) ndarray
        The matrix solution of the orthogonal Procrustes problem.

Review comment (suggestion):
I assume this? But to which function is `full_matric{e}s` a parameter?
```suggestion
    skip_backends=[("dask.array", "full_matrices=True is not supported by dask")],
```

Code after:
"""
Solve the orthogonal Procrustes problem.

"""

from scipy._lib._util import _apply_over_batch
from ._decomp_svd import svd
from scipy._lib._array_api import array_namespace, xp_capabilities, _asarray, is_numpy


__all__ = ['orthogonal_procrustes']


@xp_capabilities(
    jax_jit=False,
    skip_backends=[("dask.array", "full_matrices=True is not supported by dask")],
)
@_apply_over_batch(('A', 2), ('B', 2))
def orthogonal_procrustes(A, B, check_finite=True):
    """
    Compute the matrix solution of the orthogonal (or unitary) Procrustes problem.

    Given matrices `A` and `B` of the same shape, find an orthogonal (or unitary in
    the case of complex input) matrix `R` that most closely maps `A` to `B` using the
    algorithm given in [1]_.

    Parameters
    ----------
    A : (M, N) array_like
        Matrix to be mapped.
    B : (M, N) array_like
        Target matrix.
    check_finite : bool, optional
        Whether to check that the input matrices contain only finite numbers.
        Disabling may give a performance gain, but may result in problems
        (crashes, non-termination) if the inputs do contain infinities or NaNs.

    Returns
    -------
    R : (N, N) ndarray
        The matrix solution of the orthogonal Procrustes problem.
Repository: scipy/scipy
File: scipy/stats/_discrete_distns.py (Python)

Code before:
#
# Author:  Travis Oliphant  2002-2011 with contributions from
#          SciPy Developers 2004-2011
#
from functools import partial

from scipy import special
from scipy.special import entr, logsumexp, betaln, gammaln as gamln
from scipy.special._ufuncs import _gen_harmonic
from scipy._lib._util import rng_integers
import scipy._lib.array_api_extra as xpx
from scipy.interpolate import interp1d

from numpy import floor, ceil, log, exp, sqrt, log1p, expm1, tanh, cosh, sinh

import numpy as np

from ._distn_infrastructure import (rv_discrete, get_distribution_names,
                                    _vectorize_rvs_over_shapes,
                                    _ShapeInfo, _isintegral,
                                    rv_discrete_frozen)
from ._biasedurn import (_PyFishersNCHypergeometric,
                         _PyWalleniusNCHypergeometric,
                         _PyStochasticLib3)
from ._stats_pythran import _poisson_binom

import scipy.special._ufuncs as scu



class binom_gen(rv_discrete):
    r"""A binomial discrete random variable.

    %(before_notes)s

Review comment (suggestion):
The rest of the file seems to use imports as `scu.<fn>`. 

Code after:
#
# Author:  Travis Oliphant  2002-2011 with contributions from
#          SciPy Developers 2004-2011
#
from functools import partial

from scipy import special
from scipy.special import entr, logsumexp, betaln, gammaln as gamln
import scipy.special._ufuncs as scu
from scipy._lib._util import rng_integers
import scipy._lib.array_api_extra as xpx
from scipy.interpolate import interp1d

from numpy import floor, ceil, log, exp, sqrt, log1p, expm1, tanh, cosh, sinh

import numpy as np

from ._distn_infrastructure import (rv_discrete, get_distribution_names,
                                    _vectorize_rvs_over_shapes,
                                    _ShapeInfo, _isintegral,
                                    rv_discrete_frozen)
from ._biasedurn import (_PyFishersNCHypergeometric,
                         _PyWalleniusNCHypergeometric,
                         _PyStochasticLib3)
from ._stats_pythran import _poisson_binom


class binom_gen(rv_discrete):
    r"""A binomial discrete random variable.

    %(before_notes)s

    Notes
    -----
Repository: scipy/scipy
File: scipy/special/gen_harmonic.h (C/C++)

Code before:
        return NAN;
    }
    if (std::isinf(a)) {
        if (a > 0) {
            // a = +inf
            return 1.0;
        }
        else {
            // a = -inf
            if (n == 1) {
                // IEEE: pow(1.0, inf) = 1
                return 1.0;
            }
            else {
                return INFINITY;
            }
        }
    }
    if (a == 0) {
        return static_cast<double>(n);
    }
    if (a <= 1) {
        return sum_powers(1, n, a);
    }
    else {
        // If here, we know a is finite and a > 1, and n >= 1.
        //
        // For a > 1, we can use the formula zeta(a, 1) - zeta(a, n + 1),
        // where zeta(a, k) is the Hurwitiz zeta function.
        // But if zeta(a, 1) and zeta(a, n + 1) are close, precision is lost
        // in the subtraction, so we use the explicit sum instead. We consider
        // the values "close" if zeta(a, n + 1)/zeta(a, 1) > 0.5.
        //
        double z1 = xsf::cephes::zeta(a, static_cast<double>(1));
        double znp1 =  xsf::cephes::zeta(a, static_cast<double>(n + 1));
        if (znp1 / z1 <= 0.5) {
            return z1 - znp1;
        }
        else {
            return sum_powers(1, n, a);
        }
    }
}


//
// Computes
//
//      sum_{i=j}^{k} i**-a
//      -------------------
//      sum_{i=1}^{n} i**-a

Review comment (suggestion):
This would need to change too then.
```suggestion
        // If here, we know a is finite and a > 1, and n > 1.
```

Code after:
{
    if (n < 1) {
        sf_error("_gen_harmonic", SF_ERROR_DOMAIN,
                 "n >= 1 is required, but got n = %" PRId64, n);
        return std::numeric_limits<double>::quiet_NaN();
    }
    if (n == 1) {
        // IEEE: pow(1.0, _) is 1.0.
        return 1.0;
    }
    if (std::isnan(a)) {
        return std::numeric_limits<double>::quiet_NaN();
    }
    if (std::isinf(a)) {
        if (a > 0) {
            // a = +inf
            return 1.0;
        }
        else {
            // a = -inf
            return std::numeric_limits<double>::infinity();
        }
    }
    if (a == 0) {
        return static_cast<double>(n);
    }
    if (a <= 1) {
        return sum_powers(1, n, a);
    }
    else {
        // If here, we know a is finite and a > 1, and n > 1.
        //
        // For a > 1, we can use the formula zeta(a, 1) - zeta(a, n + 1),
        // where zeta(a, k) is the Hurwitiz zeta function.
        // But if zeta(a, 1) and zeta(a, n + 1) are close, precision is lost
        // in the subtraction, so we use the explicit sum instead. We consider
        // the values "close" if zeta(a, n + 1)/zeta(a, 1) > zeta_ratio_threshold.
        //
        double z1 = xsf::cephes::zeta(a, static_cast<double>(1));
        double znp1 =  xsf::cephes::zeta(a, static_cast<double>(n + 1));
        if (znp1 / z1 <= zeta_ratio_threshold) {
            return z1 - znp1;
        }
        else {
            return sum_powers(1, n, a);
        }
    }
}


//
Repository: scipy/scipy
File: scipy/optimize/tests/test__differential_evolution.py (Python)

Code before:
"""
Unit tests for the differential global minimization algorithm.
"""
from multiprocessing.dummy import Pool as ThreadPool
import platform
import warnings

from scipy.optimize._differentialevolution import (DifferentialEvolutionSolver,
                                                   _ConstraintWrapper)
from scipy.optimize import differential_evolution, OptimizeResult
from scipy.optimize._constraints import (Bounds, NonlinearConstraint,
                                         LinearConstraint)
from scipy.optimize import rosen, minimize
from scipy.sparse import csr_array
from scipy import stats

import numpy as np
from numpy.testing import (assert_equal, assert_allclose, assert_almost_equal,
                           assert_string_equal, assert_)
from pytest import raises as assert_raises, warns
import pytest


class TestDifferentialEvolutionSolver:

    def setup_method(self):
        self.old_seterr = np.seterr(invalid='raise')
        self.limits = np.array([[0., 0.],
                                [2., 2.]])
        self.bounds = [(0., 2.), (0., 2.)]

        self.dummy_solver = DifferentialEvolutionSolver(self.quadratic,

Review comment (bug):
`partial` is not used in this file. Maybe you added this line to the wrong file?

Code after:
"""
Unit tests for the differential global minimization algorithm.
"""
from multiprocessing.dummy import Pool as ThreadPool
import platform
import warnings
from functools import partial

from scipy.optimize._differentialevolution import (DifferentialEvolutionSolver,
                                                   _ConstraintWrapper)
from scipy.optimize import differential_evolution, OptimizeResult
from scipy.optimize._constraints import (Bounds, NonlinearConstraint,
                                         LinearConstraint)
from scipy.optimize import rosen, minimize, rosen_der
from scipy.sparse import csr_array
from scipy import stats

import numpy as np
from numpy.testing import (assert_equal, assert_allclose, assert_almost_equal,
                           assert_string_equal, assert_)
from pytest import raises as assert_raises, warns
import pytest


class TestDifferentialEvolutionSolver:

    def setup_method(self):
        self.old_seterr = np.seterr(invalid='raise')
        self.limits = np.array([[0., 0.],
                                [2., 2.]])
        self.bounds = [(0., 2.), (0., 2.)]
Repository: scipy/scipy
File: scipy/signal/_delegators.py (Python)

Code before:
    ATM, only pass tuples through. Consider updating when cupyx.lti class
    is supported.
    """
    if isinstance(arg, tuple):
        return arg
    else:
        return (None,)


def _skip_if_str_or_tuple(window):
    """Handle `window` being a str or a tuple or an array-like.
    """
    if isinstance(window, str) or isinstance(window, tuple) or callable(window):
        return None
    else:
        return window


def _skip_if_poly1d(arg):
    return None if isinstance(arg, np.poly1d) else arg


###################

def argrelextrema_signature(data, *args, **kwds):
    return array_namespace(data)

argrelmax_signature = argrelextrema_signature
argrelmin_signature = argrelextrema_signature


def band_stop_obj_signature(wp, ind, passb, stopb, gpass, gstop, type):
    return array_namespace(passb, stopb)


def bessel_signature(N, Wn, *args, **kwds):
    return array_namespace(Wn)

butter_signature = bessel_signature


def cheby2_signature(N, rs, Wn, *args, **kwds):
    return array_namespace(Wn)


def cheby1_signature(N, rp, Wn, *args, **kwds):
    return array_namespace(Wn)


def ellip_signature(N, rp, rs, Wn, *args, **kwds):
    return array_namespace(Wn)

Review comment (suggestion):
This removal makes the corresponding `delegator` to be None, meaning it will skip delegating to the `cupyx.scipy` namesake. Is this intended? I'd agree there's not much point calling into CuPy just for essentially atleast_2d.

Code after:
    ATM, only pass tuples through. Consider updating when cupyx.lti class
    is supported.
    """
    if isinstance(arg, tuple):
        return arg
    else:
        return (None,)


def _skip_if_str_or_tuple(window):
    """Handle `window` being a str or a tuple or an array-like.
    """
    if isinstance(window, str) or isinstance(window, tuple) or callable(window):
        return None
    else:
        return window


def _skip_if_poly1d(arg):
    return None if isinstance(arg, np.poly1d) else arg


###################

def abcd_normalize_signature(A=None, B=None, C=None, D=None, *, dtype=None):
    return array_namespace(A, B, C, D)

def argrelextrema_signature(data, *args, **kwds):
    return array_namespace(data)

argrelmax_signature = argrelextrema_signature
argrelmin_signature = argrelextrema_signature


def band_stop_obj_signature(wp, ind, passb, stopb, gpass, gstop, type):
    return array_namespace(passb, stopb)


def bessel_signature(N, Wn, *args, **kwds):
    return array_namespace(Wn)

butter_signature = bessel_signature


def cheby2_signature(N, rs, Wn, *args, **kwds):
    return array_namespace(Wn)


def cheby1_signature(N, rp, Wn, *args, **kwds):
    return array_namespace(Wn)
Repository: scipy/scipy
File: scipy/special/boost_special_functions.h (C/C++)

Code before:
    if ((x <= 0) || (p < 0) || (p > 1)) {
        sf_error("chdtriv", SF_ERROR_DOMAIN, NULL);
        return NAN;
    }
    if ((p == 1)) {
        return 0.0;
    }
    Real y;
    try {
        y = 2 * boost::math::gamma_p_inva<Real>(0.5 * x, p, SpecialPolicy());
    } catch (const std::domain_error& e) {
        sf_error("chdtriv", SF_ERROR_DOMAIN, NULL);
        y = NAN;
    } catch (const std::overflow_error& e) {
        sf_error("chdtriv", SF_ERROR_OVERFLOW, NULL);
        y = INFINITY;
    } catch (const std::underflow_error& e) {
        sf_error("chdtriv", SF_ERROR_UNDERFLOW, NULL);
        y = 0;
    } catch (...) {
        sf_error("chdtriv", SF_ERROR_NO_RESULT, NULL);
        y = NAN;
    }
    if (y < 0) {
        sf_error("chdtriv", SF_ERROR_NO_RESULT, NULL);
        y = NAN;
    }
    return y;
}

float
chdtriv_float(float p, float x)
{
    return chdtriv_wrap(p, x);
}

double
chdtriv_double(double p, double x)
{
    return chdtriv_wrap(p, x);
}

#endif


Review comment (suggestion):
```suggestion
    // In theory, y should never be negative, but with an overabundance of caution
    // we check and return NAN if it happens.
    if (y < 0) {
```

Code after:
    if ((x <= 0) || (p < 0) || (p > 1)) {
        sf_error("chdtriv", SF_ERROR_DOMAIN, NULL);
        return NAN;
    }
    if (p == 1) {
        return 0.0;
    }
    Real y;
    try {
        y = 2 * boost::math::gamma_p_inva<Real>(0.5 * x, p, SpecialPolicy());
    } catch (const std::domain_error& e) {
        sf_error("chdtriv", SF_ERROR_DOMAIN, NULL);
        y = NAN;
    } catch (const std::overflow_error& e) {
        sf_error("chdtriv", SF_ERROR_OVERFLOW, NULL);
        y = INFINITY;
    } catch (const std::underflow_error& e) {
        sf_error("chdtriv", SF_ERROR_UNDERFLOW, NULL);
        y = 0;
    } catch (...) {
        sf_error("chdtriv", SF_ERROR_NO_RESULT, NULL);
        y = NAN;
    }
    // In theory, y should never be negative, but with an overabundance of caution
    // we check and return NAN if it happens.
    if (y < 0) {
        sf_error("chdtriv", SF_ERROR_NO_RESULT, NULL);
        y = NAN;
    }
    return y;
}

float
chdtriv_float(float p, float x)
{
    return chdtriv_wrap(p, x);
}

double
chdtriv_double(double p, double x)
{
    return chdtriv_wrap(p, x);
}

#endif
Repository: scipy/scipy
File: doc/source/tutorial/examples/signal_Filtering_firwin2_example.py (Python)

Code before:
import matplotlib.pyplot as plt
import numpy as np
import scipy.signal as signal

fs, numtaps = 2, 41  # sampling frequency and number of taps
freqs = [0.0, 0.3, 0.6, 1.0]  # corner frequencies
gains = [1.0, 2.0, 0.5, 0.8]  # desired gains at corner frequencies

b = signal.firwin2(numtaps, freqs, gains, fs=fs)  # design filter
f, H = signal.freqz(b, fs=fs)  # calculate frequency response

fg, ax = plt.subplots(1, 1, tight_layout=True)  # do the plotting
ax.set(title=f"Frequency Response of {numtaps}-tap Band-pass FIR-Filter",
       ylabel="Gain", xlim=(0, fs/2),
       xlabel=rf"Frequency $f\,$ in hertz (sampling frequency $f_S={fs}\,$Hz)")
ax.plot(freqs,gains, 'ko--', alpha=.5, label="Desired Gains")
ax.plot(f, np.abs(H), label="Response $|H(f)|$")
ax.grid()
ax.legend()
plt.show()


Review comment (suggestion):
Optional -- but I find this to do something better 99% of the time in my own use cases nowadays
```suggestion
fg, ax = plt.subplots(1, 1, layout="constrained")  # do the plotting
```

Code after:
import matplotlib.pyplot as plt
import numpy as np
import scipy.signal as signal

fs, numtaps = 2, 41  # sampling frequency and number of taps
freqs = [0.0, 0.3, 0.6, 1.0]  # corner frequencies
gains = [1.0, 2.0, 0.5, 0.8]  # desired gains at corner frequencies

b = signal.firwin2(numtaps, freqs, gains, fs=fs)  # design filter
f, H = signal.freqz(b, fs=fs)  # calculate frequency response

fg, ax = plt.subplots(1, 1, layout="constrained")  # do the plotting
ax.set(title=f"Frequency Response of {numtaps}-tap Band-pass FIR-Filter",
       ylabel="Gain", xlim=(0, fs/2),
       xlabel=rf"Frequency $f\,$ in hertz (sampling frequency $f_S={fs}\,$Hz)")
ax.plot(freqs, gains, 'ko--', alpha=.5, label="Desired Gains")
ax.plot(f, np.abs(H), label="Response $|H(f)|$")
ax.grid()
ax.legend()
plt.show()
Repository: scipy/scipy
File: scipy/stats/_continuous_distns.py (Python)

Code before:
    See Also
    --------
    pareto

    Notes
    -----
    The probability density function for `powerlaw` is:

    .. math::

        f(x, a) = a x^{a-1}

    for :math:`0 \le x \le 1`, :math:`a > 0`.

    `powerlaw` takes ``a`` as a shape parameter for :math:`a`.

    %(after_notes)s

    For example, the support of `powerlaw` can be adjusted from the default
    interval ``[0, 1]`` to the interval ``[c, c+d]`` by setting ``loc=c`` and
    ``scale=d``. For a power-law distribution with infinite support, see
    `pareto`. For a power-law distribution described by PDF:

    .. math::

        f(x; a, l, h) = \frac{a}{h^a - l^2} x^{a-1}}

    with :math:`a \neq 0` and :math:`0 < l < x < h`, see `truncpareto`.

    `powerlaw` is a special case of `beta` with ``b=1``.

    %(example)s

    """
    def _shape_info(self):
        return [_ShapeInfo("a", False, (0, np.inf), (False, False))]

    def _pdf(self, x, a):
        # powerlaw.pdf(x, a) = a * x**(a-1)
        return a*x**(a-1.0)

    def _logpdf(self, x, a):
        return np.log(a) + sc.xlogy(a - 1, x)

    def _cdf(self, x, a):
        return x**(a*1.0)

    def _logcdf(self, x, a):
        return a*np.log(x)

    def _ppf(self, q, a):

Review comment (suggestion):
```suggestion
        f(x; a, l, h) = \frac{a}{h^a - l^2} x^{a-1}
```

Code after:

    See Also
    --------
    pareto

    Notes
    -----
    The probability density function for `powerlaw` is:

    .. math::

        f(x, a) = a x^{a-1}

    for :math:`0 \le x \le 1`, :math:`a > 0`.

    `powerlaw` takes ``a`` as a shape parameter for :math:`a`.

    %(after_notes)s

    For example, the support of `powerlaw` can be adjusted from the default
    interval ``[0, 1]`` to the interval ``[c, c+d]`` by setting ``loc=c`` and
    ``scale=d``. For a power-law distribution with infinite support, see
    `pareto`. For a power-law distribution described by PDF:

    .. math::

        f(x; a, l, h) = \frac{a}{h^a - l^2} x^{a-1}

    with :math:`a \neq 0` and :math:`0 < l < x < h`, see `truncpareto`.

    `powerlaw` is a special case of `beta` with ``b=1``.

    %(example)s

    """
    def _shape_info(self):
        return [_ShapeInfo("a", False, (0, np.inf), (False, False))]

    def _pdf(self, x, a):
        # powerlaw.pdf(x, a) = a * x**(a-1)
        return a*x**(a-1.0)

    def _logpdf(self, x, a):
        return np.log(a) + sc.xlogy(a - 1, x)

    def _cdf(self, x, a):
        return x**(a*1.0)

    def _logcdf(self, x, a):
        return a*np.log(x)
Repository: scipy/scipy
File: scipy/spatial/transform/_rotation.py (Python)

Code before:
        return out[:-1]
    return xp_promote(*args, force_floating=True, xp=xp)


class Rotation:
    """Rotation in 3 dimensions.

    This class provides an interface to initialize from and represent rotations
    with:

    - Quaternions
    - Rotation Matrices
    - Rotation Vectors
    - Modified Rodrigues Parameters
    - Euler Angles
    - Davenport Angles (Generalized Euler Angles)

    The following operations on rotations are supported:

    - Application on vectors
    - Rotation Composition
    - Rotation Inversion
    - Rotation Indexing

    A `Rotation` instance can contain a single rotation transform or rotations of
    multiple leading dimensions. E.g., it is possible to have a tensor of (N, M, K)
    rotations. When applied to other rotations or vectors, standard broadcasting rules
    apply.

    Indexing within a rotation is supported to access a subset of the rotations stored
    in a `Rotation` instance.

    To create `Rotation` objects use ``from_...`` methods (see examples below).
    ``Rotation(...)`` is not supposed to be instantiated directly.

    Attributes
    ----------
    single

    Methods
    -------
    __len__
    from_quat
    from_matrix
    from_rotvec
    from_mrp
    from_euler
    from_davenport
    as_quat
    as_matrix
    as_rotvec

Review comment (suggestion):
IMO we should call them N-dimensional arrays rather than tensors throughout, to parallel numpy arrays.

Code after:
        return out[:-1]
    return xp_promote(*args, force_floating=True, xp=xp)


class Rotation:
    """Rotation in 3 dimensions.

    This class provides an interface to initialize from and represent rotations
    with:

    - Quaternions
    - Rotation Matrices
    - Rotation Vectors
    - Modified Rodrigues Parameters
    - Euler Angles
    - Davenport Angles (Generalized Euler Angles)

    The following operations on rotations are supported:

    - Application on vectors
    - Rotation Composition
    - Rotation Inversion
    - Rotation Indexing

    A `Rotation` instance can contain a single rotation transform or rotations of
    multiple leading dimensions. E.g., it is possible to have an N-dimensional array of
    (N, M, K) rotations. When applied to other rotations or vectors, standard broadcasting
    rules apply.

    Indexing within a rotation is supported to access a subset of the rotations stored
    in a `Rotation` instance.

    To create `Rotation` objects use ``from_...`` methods (see examples below).
    ``Rotation(...)`` is not supposed to be instantiated directly.

    Attributes
    ----------
    single

    Methods
    -------
    __len__
    from_quat
    from_matrix
    from_rotvec
    from_mrp
    from_euler
    from_davenport
    as_quat
    as_matrix
    as_rotvec