Python numpy.spacing() 使用实例

The following are code examples for showing how to use . They are extracted from open source Python projects. You can vote up the examples you like or vote down the exmaples you don’t like. You can also save this page to your account.

Example 1

def test_spacing_nextafter(self):
        """Test np.spacing and np.nextafter"""
        # All non-negative finite #'s
        a = np.arange(0x7c00, dtype=uint16)
        hinf = np.array((np.inf,), dtype=float16)
        a_f16 = a.view(dtype=float16)

        assert_equal(np.spacing(a_f16[:-1]), a_f16[1:]-a_f16[:-1])

        assert_equal(np.nextafter(a_f16[:-1], hinf), a_f16[1:])
        assert_equal(np.nextafter(a_f16[0], -hinf), -a_f16[1])
        assert_equal(np.nextafter(a_f16[1:], -hinf), a_f16[:-1])

        # switch to negatives
        a |= 0x8000

        assert_equal(np.spacing(a_f16[0]), np.spacing(a_f16[1]))
        assert_equal(np.spacing(a_f16[1:]), a_f16[:-1]-a_f16[1:])

        assert_equal(np.nextafter(a_f16[0], hinf), -a_f16[1])
        assert_equal(np.nextafter(a_f16[1:], hinf), a_f16[:-1])
        assert_equal(np.nextafter(a_f16[:-1], -hinf), a_f16[1:]) 

Example 2

def validate_cov_matrix(M):
    M = (M + M.T) * 0.5
    k = 0
    I = np.eye(M.shape[0])
    while True:
        try:
            _ = np.linalg.cholesky(M)
            break
        except np.linalg.LinAlgError:
            # Find the nearest positive definite matrix for M. Modified from
            # http://www.mathworks.com/matlabcentral/fileexchange/42885-nearestspd
            # Might take several minutes
            k += 1
            w, v = np.linalg.eig(M)
            min_eig = v.min()
            M += (-min_eig * k * k + np.spacing(min_eig)) * I
    return M 

Example 3

def __init__(self, data_manager, t_layer_sizes, p_layer_sizes, dropout=0):
        print('{:25}'.format("Initializing Model"), end='', flush=True)
        self.t_layer_sizes = t_layer_sizes
        self.p_layer_sizes = p_layer_sizes
        self.dropout = dropout

        self.data_manager = data_manager
        self.t_input_size = self.data_manager.f.feature_count
        self.output_size = self.data_manager.s.information_count

        self.time_model = StackedCells(self.t_input_size, celltype=LSTM, layers = t_layer_sizes)
        self.time_model.layers.append(PassthroughLayer())

        p_input_size = t_layer_sizes[-1] + self.output_size
        self.pitch_model = StackedCells( p_input_size, celltype=LSTM, layers = p_layer_sizes)
        self.pitch_model.layers.append(Layer(p_layer_sizes[-1], self.output_size, activation = T.nnet.sigmoid))

        self.conservativity = T.fscalar()
        self.srng = T.shared_randomstreams.RandomStreams(np.random.randint(0, 1024))

        self.epsilon = np.spacing(np.float32(1.0))

        print("Done") 

Example 4

def _logL(distr, y, y_hat):
    """The log likelihood."""
    if distr in ['softplus', 'poisson']:
        eps = np.spacing(1)
        logL = np.sum(y * np.log(y_hat + eps) - y_hat)
    elif distr == 'gaussian':
        logL = -0.5 * np.sum((y - y_hat)**2)
    elif distr == 'binomial':
        # analytical formula
        logL = np.sum(y * np.log(y_hat) + (1 - y) * np.log(1 - y_hat))

        # but this prevents underflow
        # z = beta0 + np.dot(X, beta)
        # logL = np.sum(y * z - np.log(1 + np.exp(z)))
    elif distr == 'probit':
        logL = np.sum(y * np.log(y_hat) + (1 - y) * np.log(1 - y_hat))
    elif distr == 'gamma':
        # see
        # https://www.statistics.ma.tum.de/fileadmin/w00bdb/www/czado/lec8.pdf
        nu = 1.  # shape parameter, exponential for now
        logL = np.sum(nu * (-y / y_hat - np.log(y_hat)))
    return logL 

Example 5

def test_spacing_nextafter(self):
        """Test np.spacing and np.nextafter"""
        # All non-negative finite #'s
        a = np.arange(0x7c00, dtype=uint16)
        hinf = np.array((np.inf,), dtype=float16)
        a_f16 = a.view(dtype=float16)

        assert_equal(np.spacing(a_f16[:-1]), a_f16[1:]-a_f16[:-1])

        assert_equal(np.nextafter(a_f16[:-1], hinf), a_f16[1:])
        assert_equal(np.nextafter(a_f16[0], -hinf), -a_f16[1])
        assert_equal(np.nextafter(a_f16[1:], -hinf), a_f16[:-1])

        # switch to negatives
        a |= 0x8000

        assert_equal(np.spacing(a_f16[0]), np.spacing(a_f16[1]))
        assert_equal(np.spacing(a_f16[1:]), a_f16[:-1]-a_f16[1:])

        assert_equal(np.nextafter(a_f16[0], hinf), -a_f16[1])
        assert_equal(np.nextafter(a_f16[1:], hinf), a_f16[:-1])
        assert_equal(np.nextafter(a_f16[:-1], -hinf), a_f16[1:]) 

Example 6

def unitize(v):
    """
    UNIT Unitize a vector

    :param v: given unit vector
    :return: a unit-vector parallel to V.

    Reports error for the case where V is non-symbolic and norm(V) is zero
    """
    n = np.linalg.norm(v, "fro")
    # Todo ISA
    if n > np.spacing([1])[0]:
        return v / n
    else:
        raise AttributeError("Vector has zero norm")


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

Example 7

def SID(s1, s2):
    """
    Computes the spectral information divergence between two vectors.

    Parameters:
        s1: `numpy array`
            The first vector.

        s2: `numpy array`
            The second vector.

    Returns: `float`
            Spectral information divergence between s1 and s2.

    Reference
        C.-I. Chang, "An Information-Theoretic Approach to SpectralVariability,
        Similarity, and Discrimination for Hyperspectral Image"
        IEEE TRANSACTIONS ON INFORMATION THEORY, VOL. 46, NO. 5, AUGUST 2000.

    """
    p = (s1 / np.sum(s1)) + np.spacing(1)
    q = (s2 / np.sum(s2)) + np.spacing(1)
    return np.sum(p * np.log(p / q) + q * np.log(q / p)) 

Example 8

def test_spacing_nextafter(self):
        """Test np.spacing and np.nextafter"""
        # All non-negative finite #'s
        a = np.arange(0x7c00, dtype=uint16)
        hinf = np.array((np.inf,), dtype=float16)
        a_f16 = a.view(dtype=float16)

        assert_equal(np.spacing(a_f16[:-1]), a_f16[1:]-a_f16[:-1])

        assert_equal(np.nextafter(a_f16[:-1], hinf), a_f16[1:])
        assert_equal(np.nextafter(a_f16[0], -hinf), -a_f16[1])
        assert_equal(np.nextafter(a_f16[1:], -hinf), a_f16[:-1])

        # switch to negatives
        a |= 0x8000

        assert_equal(np.spacing(a_f16[0]), np.spacing(a_f16[1]))
        assert_equal(np.spacing(a_f16[1:]), a_f16[:-1]-a_f16[1:])

        assert_equal(np.nextafter(a_f16[0], hinf), -a_f16[1])
        assert_equal(np.nextafter(a_f16[1:], hinf), a_f16[:-1])
        assert_equal(np.nextafter(a_f16[:-1], -hinf), a_f16[1:]) 

Example 9

def test_spacing_nextafter(self):
        """Test np.spacing and np.nextafter"""
        # All non-negative finite #'s
        a = np.arange(0x7c00, dtype=uint16)
        hinf = np.array((np.inf,), dtype=float16)
        a_f16 = a.view(dtype=float16)

        assert_equal(np.spacing(a_f16[:-1]), a_f16[1:]-a_f16[:-1])

        assert_equal(np.nextafter(a_f16[:-1], hinf), a_f16[1:])
        assert_equal(np.nextafter(a_f16[0], -hinf), -a_f16[1])
        assert_equal(np.nextafter(a_f16[1:], -hinf), a_f16[:-1])

        # switch to negatives
        a |= 0x8000

        assert_equal(np.spacing(a_f16[0]), np.spacing(a_f16[1]))
        assert_equal(np.spacing(a_f16[1:]), a_f16[:-1]-a_f16[1:])

        assert_equal(np.nextafter(a_f16[0], hinf), -a_f16[1])
        assert_equal(np.nextafter(a_f16[1:], hinf), a_f16[:-1])
        assert_equal(np.nextafter(a_f16[:-1], -hinf), a_f16[1:]) 

Example 10

def test_spacing_nextafter(self):
        """Test np.spacing and np.nextafter"""
        # All non-negative finite #'s
        a = np.arange(0x7c00, dtype=uint16)
        hinf = np.array((np.inf,), dtype=float16)
        a_f16 = a.view(dtype=float16)

        assert_equal(np.spacing(a_f16[:-1]), a_f16[1:]-a_f16[:-1])

        assert_equal(np.nextafter(a_f16[:-1], hinf), a_f16[1:])
        assert_equal(np.nextafter(a_f16[0], -hinf), -a_f16[1])
        assert_equal(np.nextafter(a_f16[1:], -hinf), a_f16[:-1])

        # switch to negatives
        a |= 0x8000

        assert_equal(np.spacing(a_f16[0]), np.spacing(a_f16[1]))
        assert_equal(np.spacing(a_f16[1:]), a_f16[:-1]-a_f16[1:])

        assert_equal(np.nextafter(a_f16[0], hinf), -a_f16[1])
        assert_equal(np.nextafter(a_f16[1:], hinf), a_f16[:-1])
        assert_equal(np.nextafter(a_f16[:-1], -hinf), a_f16[1:]) 

Example 11

def test_spacing_nextafter(self):
        """Test np.spacing and np.nextafter"""
        # All non-negative finite #'s
        a = np.arange(0x7c00, dtype=uint16)
        hinf = np.array((np.inf,), dtype=float16)
        a_f16 = a.view(dtype=float16)

        assert_equal(np.spacing(a_f16[:-1]), a_f16[1:]-a_f16[:-1])

        assert_equal(np.nextafter(a_f16[:-1], hinf), a_f16[1:])
        assert_equal(np.nextafter(a_f16[0], -hinf), -a_f16[1])
        assert_equal(np.nextafter(a_f16[1:], -hinf), a_f16[:-1])

        # switch to negatives
        a |= 0x8000

        assert_equal(np.spacing(a_f16[0]), np.spacing(a_f16[1]))
        assert_equal(np.spacing(a_f16[1:]), a_f16[:-1]-a_f16[1:])

        assert_equal(np.nextafter(a_f16[0], hinf), -a_f16[1])
        assert_equal(np.nextafter(a_f16[1:], hinf), a_f16[:-1])
        assert_equal(np.nextafter(a_f16[:-1], -hinf), a_f16[1:]) 

Example 12

def laplacian(W, normalized=True):
    """Return the Laplacian of the weigth matrix."""

    # Degree matrix.
    d = W.sum(axis=0)

    # Laplacian matrix.
    if not normalized:
        D = scipy.sparse.diags(d.A.squeeze(), 0)
        L = D - W
    else:
        d += np.spacing(np.array(0, W.dtype))
        d = 1 / np.sqrt(d)
        D = scipy.sparse.diags(d.A.squeeze(), 0)
        I = scipy.sparse.identity(d.size, dtype=W.dtype)
        L = I - D * W * D

    # assert np.abs(L - L.T).mean() < 1e-9
    assert type(L) is scipy.sparse.csr.csr_matrix
    return L 

Example 13

def sensitivity(Ntp, Nfn, eps=numpy.spacing(1)):
    """Sensitivity

    Wikipedia entry https://en.wikipedia.org/wiki/Sensitivity_and_specificity

    Parameters
    ----------
    Ntp : int >=0
        Number of true positives

    Nfn : int >=0
        Number of false negatives

    eps : float
        eps
        (Default value=numpy.spacing(1))

    Returns
    -------
    sensitivity: float
        Sensitivity

    """

    return float(Ntp / (Ntp + Nfn + eps)) 

Example 14

def specificity(Ntn, Nfp, eps=numpy.spacing(1)):
    """Specificity

    Wikipedia entry https://en.wikipedia.org/wiki/Sensitivity_and_specificity

    Parameters
    ----------
    Ntn : int >= 0
        Number of true negatives

    Nfp : int >= 0
        Number of false positives
    
    eps : float
        eps
        (Default value=numpy.spacing(1))

    Returns
    -------
    specificity: float
        Specificity

    """

    return float(Ntn / (Ntn + Nfp + eps)) 

Example 15

def substitution_rate(Nref, Nsubstitutions, eps=numpy.spacing(1)):
    """Substitution rate

    Parameters
    ----------
    Nref : int >=0
        Number of entries in the reference

    Nsubstitutions : int >=0
        Number of substitutions

    eps : float
        eps
        (Default value=numpy.spacing(1))

    Returns
    -------
    substitution_rate: float
        Substitution rate

    """

    return float(Nsubstitutions / (Nref + eps)) 

Example 16

def deletion_rate(Nref, Ndeletions, eps=numpy.spacing(1)):
    """Deletion rate

    Parameters
    ----------
    Nref : int >=0
        Number of entries in the reference

    Ndeletions : int >=0
        Number of deletions

    eps : float
        eps
        (Default value=numpy.spacing(1))

    Returns
    -------
    deletion_rate: float
        Deletion rate
        
    """

    return float(Ndeletions / (Nref + eps)) 

Example 17

def insertion_rate(Nref, Ninsertions, eps=numpy.spacing(1)):
    """Insertion rate

    Parameters
    ----------
    Nref : int >=0
        Number of entries in the reference

    Ninsertions : int >=0
        Number of insertions

    eps : float
        eps
        (Default value=numpy.spacing(1))

    Returns
    -------
    insertion_rate: float
        Insertion rate

    """

    return float(Ninsertions / (Nref + eps)) 

Example 18

def test_spacing_nextafter(self):
        """Test np.spacing and np.nextafter"""
        # All non-negative finite #'s
        a = np.arange(0x7c00, dtype=uint16)
        hinf = np.array((np.inf,), dtype=float16)
        a_f16 = a.view(dtype=float16)

        assert_equal(np.spacing(a_f16[:-1]), a_f16[1:]-a_f16[:-1])

        assert_equal(np.nextafter(a_f16[:-1], hinf), a_f16[1:])
        assert_equal(np.nextafter(a_f16[0], -hinf), -a_f16[1])
        assert_equal(np.nextafter(a_f16[1:], -hinf), a_f16[:-1])

        # switch to negatives
        a |= 0x8000

        assert_equal(np.spacing(a_f16[0]), np.spacing(a_f16[1]))
        assert_equal(np.spacing(a_f16[1:]), a_f16[:-1]-a_f16[1:])

        assert_equal(np.nextafter(a_f16[0], hinf), -a_f16[1])
        assert_equal(np.nextafter(a_f16[1:], hinf), a_f16[:-1])
        assert_equal(np.nextafter(a_f16[:-1], -hinf), a_f16[1:]) 

Example 19

def feat_eog(signals):
    """
    calculate the EOG features
    :param signals: 1D or 2D signals
    """

    if signals.ndim == 1: signals = np.expand_dims(signals,0)
    sfreq = use_sfreq
    nsamp = float(signals.shape[1])
    w = (fft(signals,axis=1)).real   
    feats = np.zeros((signals.shape[0],15),dtype='float32')
    delta = np.sum(np.abs(w[:,np.arange(0.5*nsamp/sfreq,4*nsamp/sfreq, dtype=int)]),axis=1)
    theta = np.sum(np.abs(w[:,np.arange(4*nsamp/sfreq,8*nsamp/sfreq, dtype=int)]),axis=1)
    alpha = np.sum(np.abs(w[:,np.arange(8*nsamp/sfreq,13*nsamp/sfreq, dtype=int)]),axis=1)
    beta  = np.sum(np.abs(w[:,np.arange(13*nsamp/sfreq,20*nsamp/sfreq, dtype=int)]),axis=1)
    gamma = np.sum(np.abs(w[:,np.arange(20*nsamp/sfreq,50*nsamp/sfreq, dtype=int)]),axis=1)   # only until 50, because hz=100
    sum_abs_pow = delta + theta + alpha + beta + gamma
    feats[:,0] = delta /sum_abs_pow
    feats[:,1] = theta /sum_abs_pow
    feats[:,2] = alpha /sum_abs_pow
    feats[:,3] = beta  /sum_abs_pow
    feats[:,4] = gamma /sum_abs_pow
    feats[:,5] = np.dot(np.array([3.5,4,5,7,30]),feats[:,0:5].T ) / (sfreq/2-0.5) #smean
    feats[:,6] = np.sqrt(np.max(signals, axis=1))    #PAV
    feats[:,7] = np.sqrt(np.abs(np.min(signals, axis=1)))   #VAV   
    feats[:,8] = np.argmax(signals, axis=1)/nsamp #PAP
    feats[:,9] = np.argmin(signals, axis=1)/nsamp #VAP
    feats[:,10] = np.sqrt(np.sum(np.abs(signals), axis=1)/ np.mean(np.sum(np.abs(signals), axis=1))) # AUC
    feats[:,11] = np.sum(((np.roll(np.sign(signals), 1,axis=1) - np.sign(signals)) != 0).astype(int),axis=1)/nsamp #TVC
    feats[:,12] = np.log10(np.std(signals, axis=1)) #STD/VAR
    feats[:,13] = np.log10(stats.kurtosis(signals,fisher=False,axis=1))       # kurtosis
    feats[:,14] = np.log10(-np.sum([(x/nsamp)*((np.log((x+np.spacing(1))/nsamp))) for x in np.apply_along_axis(lambda x: np.histogram(x, bins=8)[0], 1, signals)],axis=1))  # entropy.. yay, one line...
    if np.any(feats==np.nan): print('NaN detected')
    return np.nan_to_num(feats) 

Example 20

def feat_emg(signals):
    """
    calculate the EMG median as defined by Leangkvist (2012),
    """
    if signals.ndim == 1: signals = np.expand_dims(signals,0)
    sfreq = use_sfreq
    nsamp = float(signals.shape[1])
    w = (fft(signals,axis=1)).real   
    feats = np.zeros((signals.shape[0],13),dtype='float32')
    delta = np.sum(np.abs(w[:,np.arange(0.5*nsamp/sfreq,4*nsamp/sfreq, dtype=int)]),axis=1)
    theta = np.sum(np.abs(w[:,np.arange(4*nsamp/sfreq,8*nsamp/sfreq, dtype=int)]),axis=1)
    alpha = np.sum(np.abs(w[:,np.arange(8*nsamp/sfreq,13*nsamp/sfreq, dtype=int)]),axis=1)
    beta  = np.sum(np.abs(w[:,np.arange(13*nsamp/sfreq,20*nsamp/sfreq, dtype=int)]),axis=1)
    gamma = np.sum(np.abs(w[:,np.arange(20*nsamp/sfreq,50*nsamp/sfreq, dtype=int)]),axis=1)   # only until 50, because hz=100
    sum_abs_pow = delta + theta + alpha + beta + gamma
    feats[:,0] = delta /sum_abs_pow
    feats[:,1] = theta /sum_abs_pow
    feats[:,2] = alpha /sum_abs_pow
    feats[:,3] = beta  /sum_abs_pow
    feats[:,4] = gamma /sum_abs_pow
    feats[:,5] = np.dot(np.array([3.5,4,5,7,30]),feats[:,0:5].T ) / (sfreq/2-0.5) #smean
    emg = np.sum(np.abs(w[:,np.arange(12.5*nsamp/sfreq,32*nsamp/sfreq, dtype=int)]),axis=1)
    feats[:,6] = emg / np.sum(np.abs(w[:,np.arange(8*nsamp/sfreq,32*nsamp/sfreq, dtype=int)]),axis=1)  # ratio of high freq to total motor
    feats[:,7] = np.median(np.abs(w[:,np.arange(8*nsamp/sfreq,32*nsamp/sfreq, dtype=int)]),axis=1)    # median freq
    feats[:,8] = np.mean(np.abs(w[:,np.arange(8*nsamp/sfreq,32*nsamp/sfreq, dtype=int)]),axis=1)   #  mean freq
    feats[:,9] = np.std(signals, axis=1)    #  std 
    feats[:,10] = np.mean(signals,axis=1)
    feats[:,11] = np.log10(stats.kurtosis(signals,fisher=False,axis=1) )
    feats[:,12] = np.log10(-np.sum([(x/nsamp)*((np.log((x+np.spacing(1))/nsamp))) for x in np.apply_along_axis(lambda x: np.histogram(x, bins=8)[0], 1, signals)],axis=1))  # entropy.. yay, one line...
    if np.any(feats==np.nan): print('NaN detected')

    return np.nan_to_num(feats) 

Example 21

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) 

Example 22

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_) 

Example 23

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)) 

Example 24

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()) 

Example 25

def KL(a, b):
    """Calculate the Kullback Leibler divergence between a and b """
    D_KL = np.nansum(np.multiply(a, np.log(np.divide(a, b+np.spacing(1)))), axis=1)
    return D_KL 

Example 26

def calc_information(probTgivenXs, PYgivenTs, PXs, PYs):
    """Calculate the MI - I(X;T) and I(Y;T)"""
    PTs = np.nansum(probTgivenXs*PXs, axis=1)
    Ht = np.nansum(-np.dot(PTs, np.log2(PTs)))
    Htx = - np.nansum((np.dot(np.multiply(probTgivenXs, np.log2(probTgivenXs)), PXs)))
    Hyt = - np.nansum(np.dot(PYgivenTs*np.log2(PYgivenTs+np.spacing(1)), PTs))
    Hy = np.nansum(-PYs * np.log2(PYs+np.spacing(1)))
    IYT = Hy - Hyt
    ITX = Ht - Htx
    return ITX, IYT 

Example 27

def calc_information_1(probTgivenXs, PYgivenTs, PXs, PYs, PTs):
    """Calculate the MI - I(X;T) and I(Y;T)"""
    #PTs = np.nansum(probTgivenXs*PXs, axis=1)
    Ht = np.nansum(-np.dot(PTs, np.log2(PTs+np.spacing(1))))
    Htx = - np.nansum((np.dot(np.multiply(probTgivenXs, np.log2(probTgivenXs+np.spacing(1))), PXs)))
    Hyt = - np.nansum(np.dot(PYgivenTs*np.log2(PYgivenTs+np.spacing(1)), PTs))
    Hy = np.nansum(-PYs * np.log2(PYs+np.spacing(1)))
    IYT = Hy - Hyt
    ITX = Ht - Htx
    return ITX, IYT 

Example 28

def calc_information(probTgivenXs, PYgivenTs, PXs, PYs, PTs):
    """Calculate the MI - I(X;T) and I(Y;T)"""
    #PTs = np.nansum(probTgivenXs*PXs, axis=1)
    t_indeces = np.nonzero(PTs)
    Ht = np.nansum(-np.dot(PTs, np.log2(PTs+np.spacing(1))))
    Htx = - np.nansum((np.dot(np.multiply(probTgivenXs, np.log2(probTgivenXs)), PXs)))
    Hyt = - np.nansum(np.dot(PYgivenTs*np.log2(PYgivenTs+np.spacing(1)), PTs))
    Hy = np.nansum(-PYs * np.log2(PYs+np.spacing(1)))

    IYT = Hy - Hyt
    ITX = Ht - Htx

    return ITX, IYT 

Example 29

def t_calc_information(p_x_given_t, PYgivenTs, PXs, PYs):
    """Calculate the MI - I(X;T) and I(Y;T)"""
    Hx = np.nansum(-np.dot(PXs, np.log2(PXs)))
    Hxt = - np.nansum((np.dot(np.multiply(p_x_given_t, np.log2(p_x_given_t)), PXs)))
    Hyt = - np.nansum(np.dot(PYgivenTs*np.log2(PYgivenTs+np.spacing(1)), PTs))
    Hy = np.nansum(-PYs * np.log2(PYs+np.spacing(1)))
    IYT = Hy - Hyt
    ITX = Hx - Hxt
    return ITX, IYT 

Example 30

def barycentric_coordinates_of_projection(p, q, u, v):
    """ Given a point, gives projected coords of that point to a triangle
    in barycentric coordinates.

    See Heidrich, Computing the Barycentric Coordinates of a Projected Point, JGT 05
    at http://www.cs.ubc.ca/~heidrich/Papers/JGT.05.pdf

    Args:
        p: point to project
        q: a vertex of the triangle to project into
        u,v: edges of the the triangle such that it has vertices q, q+u, q+v

    Returns:
        b: barycentric coordinates of p's projection in triangle defined by q,u,v
            vectorized so p,q,u,v can all be 3xN
    """

    p = p.T
    q = q.T
    u = u.T
    v = v.T

    n = np.cross(u, v, axis=0)
    s = np.sum(n*n, axis=0)

    # If the triangle edges are collinear, cross-product is zero,
    # which makes "s" 0, which gives us divide by zero. So we
    # make the arbitrary choice to set s to epsv (=numpy.spacing(1)),
    # the closest thing to zero
    if np.isscalar(s):
        s = s if s else np.spacing(1)
    else:
        s[s == 0] = np.spacing(1)

    oneOver4ASquared = 1.0 / s
    w = p - q
    b2 = np.sum(np.cross(u, w, axis=0) * n, axis=0) * oneOver4ASquared
    b1 = np.sum(np.cross(w, v, axis=0) * n, axis=0) * oneOver4ASquared
    b = np.vstack((1 - b1 - b2, b1, b2))

    return b.T 

Example 31

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) 

Example 32

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) 

Example 33

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_) 

Example 34

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)) 

Example 35

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()) 

Example 36

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))) 

Example 37

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])) 

Example 38

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) 

Example 39

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] 

Example 40

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


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

Example 41

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] 

Example 42

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 

Example 43

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 

Example 44

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))) 

Example 45

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) 

Example 46

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_) 

Example 47

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)) 

Example 48

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()) 

Example 49

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) 

Example 50

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 
点赞

发表评论

电子邮件地址不会被公开。 必填项已用*标注