python类spacing()的实例源码

def log_likelihood(y, yhat):
    '''Helper function to compute the log likelihood.'''
    eps = np.spacing(1)
    return np.nansum(y * np.log(eps + yhat) - yhat)

def _test_spacing(t):
    one = t(1)
    eps = np.finfo(t).eps
    nan = t(np.nan)
    inf = t(np.inf)
    with np.errstate(invalid='ignore'):
        assert_(np.spacing(one) == eps)
        assert_(np.isnan(np.spacing(nan)))
        assert_(np.isnan(np.spacing(inf)))
        assert_(np.isnan(np.spacing(-inf)))
        assert_(np.spacing(t(1e30)) != 0)

def test_spacing_gfortran():
    # Reference from this fortran file, built with gfortran 4.3.3 on linux
    # 32bits:
    #       PROGRAM test_spacing
    #        INTEGER, PARAMETER :: SGL = SELECTED_REAL_KIND(p=6, r=37)
    #        INTEGER, PARAMETER :: DBL = SELECTED_REAL_KIND(p=13, r=200)
    #
    #        WRITE(*,*) spacing(0.00001_DBL)
    #        WRITE(*,*) spacing(1.0_DBL)
    #        WRITE(*,*) spacing(1000._DBL)
    #        WRITE(*,*) spacing(10500._DBL)
    #
    #        WRITE(*,*) spacing(0.00001_SGL)
    #        WRITE(*,*) spacing(1.0_SGL)
    #        WRITE(*,*) spacing(1000._SGL)
    #        WRITE(*,*) spacing(10500._SGL)
    #       END PROGRAM
    ref = {np.float64: [1.69406589450860068E-021,
                        2.22044604925031308E-016,
                        1.13686837721616030E-013,
                        1.81898940354585648E-012],
           np.float32: [9.09494702E-13,
                        1.19209290E-07,
                        6.10351563E-05,
                        9.76562500E-04]}

    for dt, dec_ in zip([np.float32, np.float64], (10, 20)):
        x = np.array([1e-5, 1, 1000, 10500], dtype=dt)
        assert_array_almost_equal(np.spacing(x), ref[dt], decimal=dec_)

def test_nextafter_vs_spacing():
    # XXX: spacing does not handle long double yet
    for t in [np.float32, np.float64]:
        for _f in [1, 1e-5, 1000]:
            f = t(_f)
            f1 = t(_f + 1)
            assert_(np.nextafter(f, f1) - f == np.spacing(f))

def test_nans_infs(self):
        with np.errstate(all='ignore'):
            # Check some of the ufuncs
            assert_equal(np.isnan(self.all_f16), np.isnan(self.all_f32))
            assert_equal(np.isinf(self.all_f16), np.isinf(self.all_f32))
            assert_equal(np.isfinite(self.all_f16), np.isfinite(self.all_f32))
            assert_equal(np.signbit(self.all_f16), np.signbit(self.all_f32))
            assert_equal(np.spacing(float16(65504)), np.inf)

            # Check comparisons of all values with NaN
            nan = float16(np.nan)

            assert_(not (self.all_f16 == nan).any())
            assert_(not (nan == self.all_f16).any())

            assert_((self.all_f16 != nan).all())
            assert_((nan != self.all_f16).all())

            assert_(not (self.all_f16 < nan).any())
            assert_(not (nan < self.all_f16).any())

            assert_(not (self.all_f16 <= nan).any())
            assert_(not (nan <= self.all_f16).any())

            assert_(not (self.all_f16 > nan).any())
            assert_(not (nan > self.all_f16).any())

            assert_(not (self.all_f16 >= nan).any())
            assert_(not (nan >= self.all_f16).any())

def rel_entropy_true(p, q):
    """KL divergence (relative entropy) D(p||q) in bits

    Returns a scalar entropy when the input distributions p and q are
    vectors of probability masses, or returns in a row vector the
    columnwise relative entropies of the input probability matrices p and
    q"""

    if type(p) == list or type(q) == tuple:
        p = np.array(p)
    if type(q) == list or type(q) == tuple:
        q = np.array(q)

    if not p.shape == q.shape:
        raise Exception('p and q must be equal sizes',
                        'p: ' + str(p.shape),
                        'q: ' + str(q.shape))

    if (np.iscomplex(p).any() or not
        np.isfinite(p).any() or
        (p < 0).any() or
        (p > 1).any()):
        raise Exception('The probability elements of p must be real numbers'
                        'between 0 and 1.')

    if (np.iscomplex(q).any() or not
        np.isfinite(q).any() or
        (q < 0).any() or
        (q > 1).any()):
        raise Exception('The probability elements of q must be real numbers'
                        'between 0 and 1.')

    eps = math.sqrt(np.spacing(1))
    if (np.abs(np.sum(p, axis=0) - 1) > eps).any():
        raise Exception('Sum of the probability elements of p must equal 1.')
    if (np.abs(np.sum(q, axis=0) - 1) > eps).any():
        raise Exception('Sum of the probability elements of q must equal 1.')

    return sum(np.log2((p**p) / (q**p)))

def check_for_dark_states(nb, Es):
        """Check for dark states, throws a warning if it finds one."""
        dark_state_indices = np.where(np.abs(np.imag(Es)) < 10 * np.spacing(1))

        if len(dark_state_indices[0]) == 0:
            return

        import warnings
        warnings.warn('The {} block contains {} dark state(s) with generalized eigenenergie(s): {}'.format(nb, len(dark_state_indices), Es[dark_state_indices]))

def _ulps_away(value1, value2, num_bits=1):
    r"""Determines if ``value1`` is within ``n`` ULPs of ``value2``.

    Uses ``np.spacing`` to determine the unit of least precision (ULP)
    for ``value1`` and then checks that the different between the values
    does not exceed ``n`` ULPs.

    When ``value1 == 0`` or ``value2 == 0``, we instead check that the other
    is less than :math:`2^{-40}` (``_EPS``) in magnitude.

    .. note::

       There is also a Fortran implementation of this function, which
       will be used if it can be built.

    Args:
        value1 (float): The first value that being compared.
        value2 (float): The second value that being compared.
        num_bits (Optional[int]): The number of bits allowed to differ.
            Defaults to ``1``.

    Returns:
        bool: Predicate indicating if the values agree to ``n`` bits.
    """
    if value1 == 0.0:
        return abs(value2) < _EPS
    elif value2 == 0.0:
        return abs(value1) < _EPS
    else:
        local_epsilon = np.spacing(value1)  # pylint: disable=no-member
        return abs(value1 - value2) <= num_bits * abs(local_epsilon)

def two_by_two_det(mat):
    r"""Compute the determinant of a 2x2 matrix.

    .. note::

       This is used **only** by :func:`quadratic_jacobian_polynomial` and
       :func:`cubic_jacobian_polynomial`.

    This is "needed" because :func:`numpy.linalg.det` uses a more generic
    determinant implementation which can introduce rounding even when the
    simple :math:`a d - b c` will suffice. For example:

    .. doctest:: 2-by-2

       >>> import numpy as np
       >>> mat = np.asfortranarray([
       ...     [-1.5   , 0.1875],
       ...     [-1.6875, 0.0   ],
       ... ])
       >>> actual_det = -mat[0, 1] * mat[1, 0]
       >>> np_det = np.linalg.det(mat)
       >>> np.abs(actual_det - np_det) == np.spacing(actual_det)
       True

    Args:
        mat (numpy.ndarray): A 2x2 matrix.

    Returns:
        float: The determinant of ``mat``.
    """
    return mat[0, 0] * mat[1, 1] - mat[0, 1] * mat[1, 0]

def angvec2r(theta, v):
    """
    ANGVEC2R(THETA, V) is an orthonormal rotation matrix (3x3)
    equivalent to a rotation of THETA about the vector V.

    :param theta: rotation in radians
    :param v: vector
    :return: rotation matrix

    Notes::
    - If THETA == 0 then return identity matrix.
    - If THETA ~= 0 then V must have a finite length.
    """
    if np.isscalar(theta) is False or common.isvec(v) is False:
        raise AttributeError("Arguments must be theta and vector")
    # TODO implement ISA
    elif np.linalg.norm(v) < 10 * np.spacing([1])[0]:
        if False:
            raise AttributeError("Bad arguments")
        else:
            return np.eye(3)
    sk = skew(np.matrix(unitize(v)))
    m = np.eye(3) + np.sin(theta) * sk + (1 - np.cos(theta)) * sk * sk
    return m


# ---------------------------------------------------------------------------------------#

def so2_valid(obj):
    assert type(obj) is pose.SO2
    for each in obj:
        assert each.shape == (2, 2)
        assert abs(np.linalg.det(each) - 1) < np.spacing([1])[0]

def ishomog(tr, dim, rtest=''):
    """ISHOMOG Test if SE(3) homogeneous transformation matrix.
    ISHOMOG(T) is true if the argument T is of dimension 4x4 or 4x4xN, else false.
    ISHOMOG(T, 'valid') as above, but also checks the validity of the rotation sub-matrix.
    See Also: isrot, ishomog2, isvec"""

    assert dim == [3, 3] or dim == [4, 4]
    is_valid = None
    if rtest == 'valid':
        is_valid = lambda matrix: abs(np.linalg.det(matrix) - 1) < np.spacing([1])[0]
    flag = True
    if check_args.is_mat_list(tr):
        for matrix in tr:
            if not (matrix.shape[0] == dim[0] and matrix.shape[1] == dim[0]):
                flag = False
        # if rtest = 'valid'
        if flag and rtest == 'valid':
            flag = is_valid(tr[0])  # As in matlab code only first matrix is passed for validity test
            # TODO-Do we need to test all matrices in list for validity of rotation submatrix -- Yes
    elif isinstance(tr, np.matrix):
        if tr.shape[0] == dim[0] and tr.shape[1] == dim[0]:
            if flag and rtest == 'valid':
                flag = is_valid(tr)
        else:
            flag = False
    else:
        raise ValueError('Invalid data type passed to common.ishomog()')
    return flag

def SID_classifier(M, E, threshold):
    """
    Classify a HSI cube M using the spectral information divergence
    and a spectral library E.
    This function is part of the NormXCorr class.

    Parameters
        M : numpy array
          a HSI cube ((m*n) x p)

        E : numpy array
          a spectral library (N x p)

    Returns : numpy array
          a class map ((m*n))
    """
    def prob_vector_array(m):
        pv_array = np.ndarray(shape=m.shape, dtype=np.float32)
        sum_m = np.sum(m, axis=1)
        pv_array[:] = (m.T / sum_m).T
        return pv_array + np.spacing(1)

    mn = M.shape[0]
    N = E.shape[0]
    p = prob_vector_array(M)
    q = prob_vector_array(E)
    sid = np.ndarray((mn, N), dtype=np.float)
    for i in range(mn):
        pq = q[0:,:] * np.log(q[0:,:] / p[i,:])
        pp = p[i,:] * np.log(p[i,:] / q[0:,:])
        sid[i,:] = np.sum(pp[0:,:] + pq[0:,:], axis=1)
    if isinstance(threshold, float):
        cmap = _single_value_min(sid, threshold)
    elif isinstance(threshold, list):
        cmap = _multiple_values_min(sid, threshold)
    else:
        return np.argmin(sid, axis=1), sid
    return cmap, sid

def isStochastic(matrix):
    """Check that ``matrix`` is row stochastic.

    Returns
    =======
    is_stochastic : bool
        ``True`` if ``matrix`` is row stochastic, ``False`` otherwise.

    """
    try:
        absdiff = (_np.abs(matrix.sum(axis=1) - _np.ones(matrix.shape[0])))
    except AttributeError:
        matrix = _np.array(matrix)
        absdiff = (_np.abs(matrix.sum(axis=1) - _np.ones(matrix.shape[0])))
    return (absdiff.max() <= 10*_np.spacing(_np.float64(1)))

def _test_spacing(t):
    one = t(1)
    eps = np.finfo(t).eps
    nan = t(np.nan)
    inf = t(np.inf)
    with np.errstate(invalid='ignore'):
        assert_(np.spacing(one) == eps)
        assert_(np.isnan(np.spacing(nan)))
        assert_(np.isnan(np.spacing(inf)))
        assert_(np.isnan(np.spacing(-inf)))
        assert_(np.spacing(t(1e30)) != 0)

def test_spacing_gfortran():
    # Reference from this fortran file, built with gfortran 4.3.3 on linux
    # 32bits:
    #       PROGRAM test_spacing
    #        INTEGER, PARAMETER :: SGL = SELECTED_REAL_KIND(p=6, r=37)
    #        INTEGER, PARAMETER :: DBL = SELECTED_REAL_KIND(p=13, r=200)
    #
    #        WRITE(*,*) spacing(0.00001_DBL)
    #        WRITE(*,*) spacing(1.0_DBL)
    #        WRITE(*,*) spacing(1000._DBL)
    #        WRITE(*,*) spacing(10500._DBL)
    #
    #        WRITE(*,*) spacing(0.00001_SGL)
    #        WRITE(*,*) spacing(1.0_SGL)
    #        WRITE(*,*) spacing(1000._SGL)
    #        WRITE(*,*) spacing(10500._SGL)
    #       END PROGRAM
    ref = {}
    ref[np.float64] = [1.69406589450860068E-021,
           2.22044604925031308E-016,
           1.13686837721616030E-013,
           1.81898940354585648E-012]
    ref[np.float32] = [
            9.09494702E-13,
            1.19209290E-07,
            6.10351563E-05,
            9.76562500E-04]

    for dt, dec_ in zip([np.float32, np.float64], (10, 20)):
        x = np.array([1e-5, 1, 1000, 10500], dtype=dt)
        assert_array_almost_equal(np.spacing(x), ref[dt], decimal=dec_)

def test_nextafter_vs_spacing():
    # XXX: spacing does not handle long double yet
    for t in [np.float32, np.float64]:
        for _f in [1, 1e-5, 1000]:
            f = t(_f)
            f1 = t(_f + 1)
            assert_(np.nextafter(f, f1) - f == np.spacing(f))

def test_nans_infs(self):
        with np.errstate(all='ignore'):
            # Check some of the ufuncs
            assert_equal(np.isnan(self.all_f16), np.isnan(self.all_f32))
            assert_equal(np.isinf(self.all_f16), np.isinf(self.all_f32))
            assert_equal(np.isfinite(self.all_f16), np.isfinite(self.all_f32))
            assert_equal(np.signbit(self.all_f16), np.signbit(self.all_f32))
            assert_equal(np.spacing(float16(65504)), np.inf)

            # Check comparisons of all values with NaN
            nan = float16(np.nan)

            assert_(not (self.all_f16 == nan).any())
            assert_(not (nan == self.all_f16).any())

            assert_((self.all_f16 != nan).all())
            assert_((nan != self.all_f16).all())

            assert_(not (self.all_f16 < nan).any())
            assert_(not (nan < self.all_f16).any())

            assert_(not (self.all_f16 <= nan).any())
            assert_(not (nan <= self.all_f16).any())

            assert_(not (self.all_f16 > nan).any())
            assert_(not (nan > self.all_f16).any())

            assert_(not (self.all_f16 >= nan).any())
            assert_(not (nan >= self.all_f16).any())

def pinv(x):
    #return np.linalg.pinv(x, max(x.shape) * np.spacing(np.linalg.norm(x)))
    #return scipy.linalg.pinv(x, max(x.shape) * np.spacing(np.linalg.norm(x)))
    return scipy.linalg.pinv(x, 1e-6)

def weights(self, nlags):
        """ Evenly-spaced beta weights
        """
        eps = np.spacing(1)
        u = np.linspace(eps, 1.0 - eps, nlags)

        beta_vals = u ** (self.theta1 - 1) * (1 - u) ** (self.theta2 - 1)

        beta_vals = beta_vals / sum(beta_vals)

        if self.theta3 is not None:
            w = beta_vals + self.theta3
            return w / sum(w)

        return beta_vals