Python scipy.special.poch() Examples

def torontonian_analytical(l, nbar):
    r"""Return the value of the Torontonian of the O matrices generated by gen_omats

    Args:
        l (int): number of modes
        nbar (float): mean photon number of the first mode (the only one not prepared in vacuum)

    Returns:
        float: Value of the torontonian of gen_omats(l,nbar)
    """
    if np.allclose(l, nbar, atol=1e-14, rtol=0.0):
        return 1.0
    beta = -(nbar / (l * (1 + nbar)))
    pref = factorial(l) / beta
    p1 = pref * l / poch(1 / beta, l + 2)
    p2 = pref * beta / poch(2 + 1 / beta, l)
    return (p1 + p2) * (-1) ** l 

def test_rf(self):
        if LooseVersion(mpmath.__version__) >= LooseVersion("1.0.0"):
            # no workarounds needed
            mppoch = mpmath.rf
        else:
            def mppoch(a, m):
                # deal with cases where the result in double precision
                # hits exactly a non-positive integer, but the
                # corresponding extended-precision mpf floats don't
                if float(a + m) == int(a + m) and float(a + m) <= 0:
                    a = mpmath.mpf(a)
                    m = int(a + m) - a
                return mpmath.rf(a, m)

        assert_mpmath_equal(sc.poch,
                            mppoch,
                            [Arg(), Arg()],
                            dps=400) 

def _derivative_inplace(self, nu, axis):
        """
        Compute 1D derivative along a selected dimension in-place
        May result to non-contiguous c array.
        """
        if nu < 0:
            return self._antiderivative_inplace(-nu, axis)

        ndim = len(self.x)
        axis = axis % ndim

        # reduce order
        if nu == 0:
            # noop
            return
        else:
            sl = [slice(None)]*ndim
            sl[axis] = slice(None, -nu, None)
            c2 = self.c[sl]

        if c2.shape[axis] == 0:
            # derivative of order 0 is zero
            shp = list(c2.shape)
            shp[axis] = 1
            c2 = np.zeros(shp, dtype=c2.dtype)

        # multiply by the correct rising factorials
        factor = spec.poch(np.arange(c2.shape[axis], 0, -1), nu)
        sl = [None]*c2.ndim
        sl[axis] = slice(None)
        c2 *= factor[sl]

        self.c = c2 

def _munp(self, n, a, c):
        # Pochhammer symbol: sc.pocha,n) = gamma(a+n)/gamma(a)
        return sc.poch(a, n*1.0/c) 

def _antiderivative_inplace(self, nu, axis):
        """
        Compute 1D antiderivative along a selected dimension
        May result to non-contiguous c array.
        """
        if nu <= 0:
            return self._derivative_inplace(-nu, axis)

        ndim = len(self.x)
        axis = axis % ndim

        perm = list(range(ndim))
        perm[0], perm[axis] = perm[axis], perm[0]
        perm = perm + list(range(ndim, self.c.ndim))

        c = self.c.transpose(perm)

        c2 = np.zeros((c.shape[0] + nu,) + c.shape[1:],
                     dtype=c.dtype)
        c2[:-nu] = c

        # divide by the correct rising factorials
        factor = spec.poch(np.arange(c.shape[0], 0, -1), nu)
        c2[:-nu] /= factor[(slice(None),) + (None,)*(c.ndim-1)]

        # fix continuity of added degrees of freedom
        perm2 = list(range(c2.ndim))
        perm2[1], perm2[ndim+axis] = perm2[ndim+axis], perm2[1]

        c2 = c2.transpose(perm2)
        c2 = c2.copy()
        _ppoly.fix_continuity(c2.reshape(c2.shape[0], c2.shape[1], -1),
                              self.x[axis], nu-1)

        c2 = c2.transpose(perm2)
        c2 = c2.transpose(perm)

        # Done
        self.c = c2 

def _derivative_inplace(self, nu, axis):
        """
        Compute 1D derivative along a selected dimension in-place
        May result to non-contiguous c array.
        """
        if nu < 0:
            return self._antiderivative_inplace(-nu, axis)

        ndim = len(self.x)
        axis = axis % ndim

        # reduce order
        if nu == 0:
            # noop
            return
        else:
            sl = [slice(None)]*ndim
            sl[axis] = slice(None, -nu, None)
            c2 = self.c[sl]

        if c2.shape[axis] == 0:
            # derivative of order 0 is zero
            shp = list(c2.shape)
            shp[axis] = 1
            c2 = np.zeros(shp, dtype=c2.dtype)

        # multiply by the correct rising factorials
        factor = spec.poch(np.arange(c2.shape[axis], 0, -1), nu)
        sl = [None]*c2.ndim
        sl[axis] = slice(None)
        c2 *= factor[sl]

        self.c = c2 

def _antiderivative_inplace(self, nu, axis):
        """
        Compute 1D antiderivative along a selected dimension
        May result to non-contiguous c array.
        """
        if nu <= 0:
            return self._derivative_inplace(-nu, axis)

        ndim = len(self.x)
        axis = axis % ndim

        perm = list(range(ndim))
        perm[0], perm[axis] = perm[axis], perm[0]
        perm = perm + list(range(ndim, self.c.ndim))

        c = self.c.transpose(perm)

        c2 = np.zeros((c.shape[0] + nu,) + c.shape[1:],
                     dtype=c.dtype)
        c2[:-nu] = c

        # divide by the correct rising factorials
        factor = spec.poch(np.arange(c.shape[0], 0, -1), nu)
        c2[:-nu] /= factor[(slice(None),) + (None,)*(c.ndim-1)]

        # fix continuity of added degrees of freedom
        perm2 = list(range(c2.ndim))
        perm2[1], perm2[ndim+axis] = perm2[ndim+axis], perm2[1]

        c2 = c2.transpose(perm2)
        c2 = c2.copy()
        _ppoly.fix_continuity(c2.reshape(c2.shape[0], c2.shape[1], -1),
                              self.x[axis], nu-1)

        c2 = c2.transpose(perm2)
        c2 = c2.transpose(perm)

        # Done
        self.c = c2 

def _dpow(x, y, n):
    """
    d^n (x**y) / dx^n
    """
    if n < 0:
        raise ValueError("invalid derivative order")
    elif n > y:
        return 0
    else:
        return poch(y - n + 1, n) * x**(y - n) 

def test_rf(self):
        assert_mpmath_equal(sc.poch,
                            mpmath.rf,
                            [Arg(), Arg()], dps=100) 

def _munp(self, n, a, c):
        # Pochhammer symbol: sc.pocha,n) = gamma(a+n)/gamma(a)
        return sc.poch(a, n*1.0/c) 

def _antiderivative_inplace(self, nu, axis):
        """
        Compute 1D antiderivative along a selected dimension
        May result to non-contiguous c array.
        """
        if nu <= 0:
            return self._derivative_inplace(-nu, axis)

        ndim = len(self.x)
        axis = axis % ndim

        perm = list(range(ndim))
        perm[0], perm[axis] = perm[axis], perm[0]
        perm = perm + list(range(ndim, self.c.ndim))

        c = self.c.transpose(perm)

        c2 = np.zeros((c.shape[0] + nu,) + c.shape[1:],
                     dtype=c.dtype)
        c2[:-nu] = c

        # divide by the correct rising factorials
        factor = spec.poch(np.arange(c.shape[0], 0, -1), nu)
        c2[:-nu] /= factor[(slice(None),) + (None,)*(c.ndim-1)]

        # fix continuity of added degrees of freedom
        perm2 = list(range(c2.ndim))
        perm2[1], perm2[ndim+axis] = perm2[ndim+axis], perm2[1]

        c2 = c2.transpose(perm2)
        c2 = c2.copy()
        _ppoly.fix_continuity(c2.reshape(c2.shape[0], c2.shape[1], -1),
                              self.x[axis], nu-1)

        c2 = c2.transpose(perm2)
        c2 = c2.transpose(perm)

        # Done
        self.c = c2 

def _derivative_inplace(self, nu, axis):
        """
        Compute 1D derivative along a selected dimension in-place
        May result to non-contiguous c array.
        """
        if nu < 0:
            return self._antiderivative_inplace(-nu, axis)

        ndim = len(self.x)
        axis = axis % ndim

        # reduce order
        if nu == 0:
            # noop
            return
        else:
            sl = [slice(None)]*ndim
            sl[axis] = slice(None, -nu, None)
            c2 = self.c[sl]

        if c2.shape[axis] == 0:
            # derivative of order 0 is zero
            shp = list(c2.shape)
            shp[axis] = 1
            c2 = np.zeros(shp, dtype=c2.dtype)

        # multiply by the correct rising factorials
        factor = spec.poch(np.arange(c2.shape[axis], 0, -1), nu)
        sl = [None]*c2.ndim
        sl[axis] = slice(None)
        c2 *= factor[sl]

        self.c = c2 

def derivative(self, nu=1):
        """
        Construct a new piecewise polynomial representing the derivative.

        Parameters
        ----------
        nu : int, optional
            Order of derivative to evaluate. Default is 1, i.e. compute the
            first derivative. If negative, the antiderivative is returned.

        Returns
        -------
        pp : PPoly
            Piecewise polynomial of order k2 = k - n representing the derivative
            of this polynomial.

        Notes
        -----
        Derivatives are evaluated piecewise for each polynomial
        segment, even if the polynomial is not differentiable at the
        breakpoints. The polynomial intervals are considered half-open,
        ``[a, b)``, except for the last interval which is closed
        ``[a, b]``.
        """
        if nu < 0:
            return self.antiderivative(-nu)

        # reduce order
        if nu == 0:
            c2 = self.c.copy()
        else:
            c2 = self.c[:-nu, :].copy()

        if c2.shape[0] == 0:
            # derivative of order 0 is zero
            c2 = np.zeros((1,) + c2.shape[1:], dtype=c2.dtype)

        # multiply by the correct rising factorials
        factor = spec.poch(np.arange(c2.shape[0], 0, -1), nu)
        c2 *= factor[(slice(None),) + (None,)*(c2.ndim-1)]

        # construct a compatible polynomial
        return self.construct_fast(c2, self.x, self.extrapolate, self.axis) 

def antiderivative(self, nu=1):
        """
        Construct a new piecewise polynomial representing the antiderivative.

        Antiderivative is also the indefinite integral of the function,
        and derivative is its inverse operation.

        Parameters
        ----------
        nu : int, optional
            Order of antiderivative to evaluate. Default is 1, i.e. compute
            the first integral. If negative, the derivative is returned.

        Returns
        -------
        pp : PPoly
            Piecewise polynomial of order k2 = k + n representing
            the antiderivative of this polynomial.

        Notes
        -----
        The antiderivative returned by this function is continuous and
        continuously differentiable to order n-1, up to floating point
        rounding error.

        If antiderivative is computed and ``self.extrapolate='periodic'``,
        it will be set to False for the returned instance. This is done because
        the antiderivative is no longer periodic and its correct evaluation
        outside of the initially given x interval is difficult.
        """
        if nu <= 0:
            return self.derivative(-nu)

        c = np.zeros((self.c.shape[0] + nu, self.c.shape[1]) + self.c.shape[2:],
                     dtype=self.c.dtype)
        c[:-nu] = self.c

        # divide by the correct rising factorials
        factor = spec.poch(np.arange(self.c.shape[0], 0, -1), nu)
        c[:-nu] /= factor[(slice(None),) + (None,)*(c.ndim-1)]

        # fix continuity of added degrees of freedom
        self._ensure_c_contiguous()
        _ppoly.fix_continuity(c.reshape(c.shape[0], c.shape[1], -1),
                              self.x, nu - 1)

        if self.extrapolate == 'periodic':
            extrapolate = False
        else:
            extrapolate = self.extrapolate

        # construct a compatible polynomial
        return self.construct_fast(c, self.x, extrapolate, self.axis) 

def derivative(self, nu=1):
        """
        Construct a new piecewise polynomial representing the derivative.

        Parameters
        ----------
        nu : int, optional
            Order of derivative to evaluate. Default is 1, i.e. compute the
            first derivative. If negative, the antiderivative is returned.

        Returns
        -------
        pp : PPoly
            Piecewise polynomial of order k2 = k - n representing the derivative
            of this polynomial.

        Notes
        -----
        Derivatives are evaluated piecewise for each polynomial
        segment, even if the polynomial is not differentiable at the
        breakpoints. The polynomial intervals are considered half-open,
        ``[a, b)``, except for the last interval which is closed
        ``[a, b]``.
        """
        if nu < 0:
            return self.antiderivative(-nu)

        # reduce order
        if nu == 0:
            c2 = self.c.copy()
        else:
            c2 = self.c[:-nu, :].copy()

        if c2.shape[0] == 0:
            # derivative of order 0 is zero
            c2 = np.zeros((1,) + c2.shape[1:], dtype=c2.dtype)

        # multiply by the correct rising factorials
        factor = spec.poch(np.arange(c2.shape[0], 0, -1), nu)
        c2 *= factor[(slice(None),) + (None,)*(c2.ndim-1)]

        # construct a compatible polynomial
        return self.construct_fast(c2, self.x, self.extrapolate, self.axis) 

def derivative(self, nu=1):
        """
        Construct a new piecewise polynomial representing the derivative.

        Parameters
        ----------
        nu : int, optional
            Order of derivative to evaluate. Default is 1, i.e. compute the
            first derivative. If negative, the antiderivative is returned.

        Returns
        -------
        pp : PPoly
            Piecewise polynomial of order k2 = k - n representing the derivative
            of this polynomial.

        Notes
        -----
        Derivatives are evaluated piecewise for each polynomial
        segment, even if the polynomial is not differentiable at the
        breakpoints. The polynomial intervals are considered half-open,
        ``[a, b)``, except for the last interval which is closed
        ``[a, b]``.
        """
        if nu < 0:
            return self.antiderivative(-nu)

        # reduce order
        if nu == 0:
            c2 = self.c.copy()
        else:
            c2 = self.c[:-nu, :].copy()

        if c2.shape[0] == 0:
            # derivative of order 0 is zero
            c2 = np.zeros((1,) + c2.shape[1:], dtype=c2.dtype)

        # multiply by the correct rising factorials
        factor = spec.poch(np.arange(c2.shape[0], 0, -1), nu)
        c2 *= factor[(slice(None),) + (None,)*(c2.ndim-1)]

        # construct a compatible polynomial
        return self.construct_fast(c2, self.x, self.extrapolate, self.axis) 

def antiderivative(self, nu=1):
        """
        Construct a new piecewise polynomial representing the antiderivative.

        Antiderivative is also the indefinite integral of the function,
        and derivative is its inverse operation.

        Parameters
        ----------
        nu : int, optional
            Order of antiderivative to evaluate. Default is 1, i.e. compute
            the first integral. If negative, the derivative is returned.

        Returns
        -------
        pp : PPoly
            Piecewise polynomial of order k2 = k + n representing
            the antiderivative of this polynomial.

        Notes
        -----
        The antiderivative returned by this function is continuous and
        continuously differentiable to order n-1, up to floating point
        rounding error.

        If antiderivative is computed and ``self.extrapolate='periodic'``,
        it will be set to False for the returned instance. This is done because
        the antiderivative is no longer periodic and its correct evaluation
        outside of the initially given x interval is difficult.
        """
        if nu <= 0:
            return self.derivative(-nu)

        c = np.zeros((self.c.shape[0] + nu, self.c.shape[1]) + self.c.shape[2:],
                     dtype=self.c.dtype)
        c[:-nu] = self.c

        # divide by the correct rising factorials
        factor = spec.poch(np.arange(self.c.shape[0], 0, -1), nu)
        c[:-nu] /= factor[(slice(None),) + (None,)*(c.ndim-1)]

        # fix continuity of added degrees of freedom
        self._ensure_c_contiguous()
        _ppoly.fix_continuity(c.reshape(c.shape[0], c.shape[1], -1),
                              self.x, nu - 1)

        if self.extrapolate == 'periodic':
            extrapolate = False
        else:
            extrapolate = self.extrapolate

        # construct a compatible polynomial
        return self.construct_fast(c, self.x, extrapolate, self.axis) 

def antiderivative(self, nu=1):
        """
        Construct a new piecewise polynomial representing the antiderivative.

        Antiderivative is also the indefinite integral of the function,
        and derivative is its inverse operation.

        Parameters
        ----------
        nu : int, optional
            Order of antiderivative to evaluate. Default is 1, i.e. compute
            the first integral. If negative, the derivative is returned.

        Returns
        -------
        pp : PPoly
            Piecewise polynomial of order k2 = k + n representing
            the antiderivative of this polynomial.

        Notes
        -----
        The antiderivative returned by this function is continuous and
        continuously differentiable to order n-1, up to floating point
        rounding error.

        If antiderivative is computed and ``self.extrapolate='periodic'``,
        it will be set to False for the returned instance. This is done because
        the antiderivative is no longer periodic and its correct evaluation
        outside of the initially given x interval is difficult.
        """
        if nu <= 0:
            return self.derivative(-nu)

        c = np.zeros((self.c.shape[0] + nu, self.c.shape[1]) + self.c.shape[2:],
                     dtype=self.c.dtype)
        c[:-nu] = self.c

        # divide by the correct rising factorials
        factor = spec.poch(np.arange(self.c.shape[0], 0, -1), nu)
        c[:-nu] /= factor[(slice(None),) + (None,)*(c.ndim-1)]

        # fix continuity of added degrees of freedom
        self._ensure_c_contiguous()
        _ppoly.fix_continuity(c.reshape(c.shape[0], c.shape[1], -1),
                              self.x, nu - 1)

        if self.extrapolate == 'periodic':
            extrapolate = False
        else:
            extrapolate = self.extrapolate

        # construct a compatible polynomial
        return self.construct_fast(c, self.x, extrapolate, self.axis)