Python numpy.finfo函数代码示例

    def compute_nmi(self, z1, z2):
        # compute normalized mutual information
        
        n = np.size(z1)
        k1 = np.size(np.unique(z1))
        k2 = np.size(np.unique(z2))
        
        nk1 = np.zeros((k1,1))
        nk2 = np.zeros((k2,1))

        for kk in range(k1):
            nk1[kk] = np.sum(z1==kk)
        for kk in range(k2):
            nk2[kk] = np.sum(z2==kk)
            
        pk1 = nk1/float(np.sum(nk1))
        pk2 = nk2/float(np.sum(nk2))
        
        nk12 = np.zeros((k1,k2))
        for ii in range(k1):
            for jj in range(k2):
                nk12[ii,jj] = np.sum((z1==ii)*(z2==jj))
        pk12 = nk12/float(n)        
        
        Hx = -np.sum(pk1 * np.log(pk1 + np.finfo(float).eps))
        Hy = -np.sum(pk2 * np.log(pk2 + np.finfo(float).eps))
        
        Hxy = -np.sum(pk12 * np.log(pk12 + np.finfo(float).eps))
        
        MI = Hx + Hy - Hxy;
        nmi = MI/float(0.5*(Hx+Hy))
        
        return nmi

def recognition(obs, oracle, order=1, smooth=False):
    hmm_tensor = extract_hmm_tensor(oracle, max_lrs=order, smooth=smooth)

    cluster_means = np.array([np.median(oracle.f_array.data[np.array(c), :].T, axis=1)
                              for c in oracle.latent])

    cluster_means += np.finfo('float').eps
    cluster_means = (cluster_means.T / np.sum(cluster_means, axis=1)).T

    a = hmm_tensor[-1]
    a += np.finfo('float').eps
    a += 1.0
    divider = np.sum(a, axis=1)
    a = np.divide(a.T, divider).T
    log_a = np.log(a)
    # hist = np.array([len(c) for c in oracle.latent])/float(oracle.n_states-1)

    v = np.zeros((len(obs), len(oracle.latent)))
    p = np.zeros(v.shape)
    v[0] = np.log(np.dot(cluster_means, obs[0])) + np.log(1.0/len(oracle.latent))
    # v[0] = np.log(np.dot(cluster_means, obs[0])) + np.log(hist)
    # p[0] = np.arange(len(oracle.latent))
    for t in range(1, len(obs)):
        s = v[t-1]+log_a.T
        v[t] = np.max(s, axis=1)+np.log(np.dot(cluster_means, obs[t]))
        p[t-1] = np.argmax(s, axis=1)

    return v, p

def _gpinv(p, k, sigma):
    """Inverse Generalized Pareto distribution function"""
    x = np.full_like(p, np.nan)
    if sigma <= 0:
        return x
    ok = (p > 0) & (p < 1)
    if np.all(ok):
        if np.abs(k) < np.finfo(float).eps:
            x = - np.log1p(-p)
        else:
            x = np.expm1(-k * np.log1p(-p)) / k
        x *= sigma
    else:
        if np.abs(k) < np.finfo(float).eps:
            x[ok] = - np.log1p(-p[ok])
        else:
            x[ok] = np.expm1(-k * np.log1p(-p[ok])) / k
        x *= sigma
        x[p == 0] = 0
        if k >= 0:
            x[p == 1] = np.inf
        else:
            x[p == 1] = - sigma / k

    return x

def quatAverage(q_in, qsym):
    """
    """
    assert q_in.ndim == 2, 'input must be 2-s hstacked quats'

    # renormalize
    q_in = unitVector(q_in)

    # check to see num of quats is > 1
    if q_in.shape[1] < 3:
        if q_in.shape[1] == 1:
            q_bar = q_in
        else:
            ma, mq = misorientation(q_in[:, 0].reshape(4, 1),
                                    q_in[:, 1].reshape(4, 1), (qsym,))
            q_bar = quatProduct(q_in[:, 0].reshape(4, 1),
                                quatOfExpMap(0.5*ma*unitVector(mq[1:].reshape(3, 1))))
    else:
        # use first quat as initial guess
        phi = 2. * arccos(q_in[0, 0])
        if phi <= finfo(float).eps:
            x0 = zeros(3)
        else:
            n = unitVector(q_in[1:, 0].reshape(3, 1))
            x0 = phi*n.flatten()
        results = leastsq(quatAverage_obj, x0, args=(q_in, qsym))
        phi = sqrt(sum(results[0]*results[0]))
        if phi <= finfo(float).eps:
            q_bar = c_[1., 0., 0., 0.].T
        else:
            n     = results[0] / phi
            q_bar = hstack([cos(0.5*phi), sin(0.5*phi)*n]).reshape(4, 1)
    return q_bar

 def test_nan(self):
     # Test that nan is 'far' from small, tiny, inf, max and min
     for dt in [np.float32, np.float64]:
         if dt == np.float32:
             maxulp = 1e6
         else:
             maxulp = 1e12
         inf = np.array([np.inf]).astype(dt)
         nan = np.array([np.nan]).astype(dt)
         big = np.array([np.finfo(dt).max])
         tiny = np.array([np.finfo(dt).tiny])
         zero = np.array([np.PZERO]).astype(dt)
         nzero = np.array([np.NZERO]).astype(dt)
         self.assertRaises(AssertionError,
                               lambda: assert_array_max_ulp(nan, inf,
                                                            maxulp=maxulp))
         self.assertRaises(AssertionError,
                               lambda: assert_array_max_ulp(nan, big,
                                                            maxulp=maxulp))
         self.assertRaises(AssertionError,
                               lambda: assert_array_max_ulp(nan, tiny,
                                                            maxulp=maxulp))
         self.assertRaises(AssertionError,
                               lambda: assert_array_max_ulp(nan, zero,
                                                            maxulp=maxulp))
         self.assertRaises(AssertionError,
                               lambda: assert_array_max_ulp(nan, nzero,
                                                            maxulp=maxulp))

def test_ulp():
    assert_equal(ulp(), np.finfo(np.float64).eps)
    assert_equal(ulp(1.0), np.finfo(np.float64).eps)
    assert_equal(ulp(np.float32(1.0)), np.finfo(np.float32).eps)
    assert_equal(ulp(np.float32(1.999)), np.finfo(np.float32).eps)
    # Integers always return 1
    assert_equal(ulp(1), 1)
    assert_equal(ulp(2**63-1), 1)
    # negative / positive same
    assert_equal(ulp(-1), 1)
    assert_equal(ulp(7.999), ulp(4.0))
    assert_equal(ulp(-7.999), ulp(4.0))
    assert_equal(ulp(np.float64(2**54-2)), 2)
    assert_equal(ulp(np.float64(2**54)), 4)
    assert_equal(ulp(np.float64(2**54)), 4)
    # Infs, NaNs return nan
    assert_true(np.isnan(ulp(np.inf)))
    assert_true(np.isnan(ulp(-np.inf)))
    assert_true(np.isnan(ulp(np.nan)))
    # 0 gives subnormal smallest
    subn64 = np.float64(2**(-1022-52))
    subn32 = np.float32(2**(-126-23))
    assert_equal(ulp(0.0), subn64)
    assert_equal(ulp(np.float64(0)), subn64)
    assert_equal(ulp(np.float32(0)), subn32)
    # as do multiples of subnormal smallest
    assert_equal(ulp(subn64 * np.float64(2**52)), subn64)
    assert_equal(ulp(subn64 * np.float64(2**53)), subn64*2)
    assert_equal(ulp(subn32 * np.float32(2**23)), subn32)
    assert_equal(ulp(subn32 * np.float32(2**24)), subn32*2)

    def _hessian_main(self, params):
        params_infl = params[:self.k_inflate]
        params_main = params[self.k_inflate:]

        y = self.endog
        w = self.model_infl.predict(params_infl)
        w = np.clip(w, np.finfo(float).eps, 1 - np.finfo(float).eps)
        score = self.score(params)
        zero_idx = np.nonzero(y == 0)[0]
        nonzero_idx = np.nonzero(y)[0]

        mu = self.model_main.predict(params_main)

        hess_arr = np.zeros((self.k_exog, self.k_exog))

        coeff = (1 + w[zero_idx] * (np.exp(mu[zero_idx]) - 1))

        #d2l/dp2
        for i in range(self.k_exog):
            for j in range(i, -1, -1):
                hess_arr[i, j] = ((
                    self.exog[zero_idx, i] * self.exog[zero_idx, j] *
                    mu[zero_idx] * (w[zero_idx] - 1) * (1 / coeff -
                    w[zero_idx] * mu[zero_idx] * np.exp(mu[zero_idx]) /
                    coeff**2)).sum() - (mu[nonzero_idx] * self.exog[nonzero_idx, i] *
                    self.exog[nonzero_idx, j]).sum())

        return hess_arr

def system_diagnostic():
    """ return various helpful/informative information about the
    current system. For instance versions of python & available packages.

    Mostly undocumented function...
    """

    # There is probably a more clever way to do the following via introspection?

    import platform
    import os
    import poppy
    import numpy
    from .version import version
    try:
        import ttk
        ttk_version = ttk.__version__
    except ImportError:
        ttk_version = 'not found'

    try:
        import wx
        wx_version = wx.__version__
    except ImportError:
        wx_version = 'not found'

    try:
        import pyfftw
        pyfftw_version = pyfftw.version
    except ImportError:
        pyfftw_version = 'not found'

    try:
        import pysynphot
        pysynphot_version = pysynphot.__version__
    except ImportError:
        pysynphot_version = 'not found'


    try:
        import astropy
        astropy_version = astropy.__version__
    except ImportError:
        astropy_version = 'not found'

    result = DIAGNOSTIC_REPORT.format(
        os=platform.platform(),
        numpy=numpy.__version__,
        python=sys.version.replace("\n", " "),
        poppy=poppy.__version__,
        webbpsf=version,
        tkinter=ttk_version,
        wxpython=wx_version,
        pyfftw=pyfftw_version,
        pysyn=pysynphot_version,
        astropy=astropy_version,
        finfo_float=numpy.finfo(numpy.float),
        finfo_complex=numpy.finfo(numpy.complex),
    )
    return result

def get_matrix_2d_ragged(workspace, distribution, histogram2D=False, transpose=False):
    num_hist = workspace.getNumberHistograms()
    delta = numpy.finfo(numpy.float64).max
    min_value = numpy.finfo(numpy.float64).max
    max_value = numpy.finfo(numpy.float64).min
    for i in range(num_hist):
        xtmp = workspace.readX(i)
        if workspace.isHistogramData():
            #input x is edges
            xtmp = mantid.plots.helperfunctions.points_from_boundaries(xtmp)
        else:
            #input x is centers
            pass
        min_value = min(min_value, xtmp.min())
        max_value = max(max_value, xtmp.max())
        diff = xtmp[1:] - xtmp[:-1]
        delta = min(delta, diff.min())
    num_edges = int(numpy.ceil((max_value - min_value)/delta)) + 1
    x_centers = numpy.linspace(min_value, max_value, num=num_edges)
    y = mantid.plots.helperfunctions.boundaries_from_points(workspace.getAxis(1).extractValues())
    z = numpy.empty([num_hist, num_edges], dtype=numpy.float64)
    for i in range(num_hist):
        centers, ztmp, _, _ = mantid.plots.helperfunctions.get_spectrum(workspace, i, distribution=distribution, withDy=False, withDx=False)
        f = interp1d(centers, ztmp, kind='nearest', bounds_error=False, fill_value=numpy.nan)
        z[i] = f(x_centers)
    if histogram2D:
        x = mantid.plots.helperfunctions.boundaries_from_points(x_centers)
    else:
        x = x_centers
    if transpose:
        return y.T,x.T,z.T
    else:
        return x,y,z

    def __init__(self):
        self.gravity = 9.8
        self.masscart = 1.0
        self.masspole = 0.1
        self.total_mass = (self.masspole + self.masscart)
        self.length = 0.5 # actually half the pole's length
        self.polemass_length = (self.masspole * self.length)
        self.force_mag = 10.0
        self.tau = 0.02  # seconds between state updates

        # Angle at which to fail the episode
        self.theta_threshold_radians = 12 * 2 * math.pi / 360
        self.x_threshold = 2.4

        # Angle limit set to 2 * theta_threshold_radians so failing observation is still within bounds
        high = np.array([
            self.x_threshold * 2,
            np.finfo(np.float32).max,
            self.theta_threshold_radians * 2,
            np.finfo(np.float32).max])

        self.action_space = spaces.Discrete(2)
        self.observation_space = spaces.Box(-high, high)

        self.seed()
        self.viewer = None
        self.state = None

        self.steps_beyond_done = None

 def almost_equal(self, other,
                  rtol=2**4 * np.finfo(np.zeros(1).dtype).eps,
                  atol=0.):
     # assert atol == 0., 'Not supported'
     assert type(other) == type(self)
     return self._impl.almost_equal(other._impl,
                                    rtol if rtol is not None else 2**4 *np.finfo(np.zeros(1.).dtype).eps)

    def get_cost(self, weights, data, display=False):
        # print 'getting cost...'
        N = float(len(data))
        reg = (self.lmbda / (2.0 * N)) * np.sum(weights ** 2)
        # reg = (self.lmbda / self.N) * np.sum(np.abs(weights))
        # self.set_network_weights(weights)
        layers = self.convert_weights_to_layers(weights)
        cost = 0.0
        for d, l in data[:]:
            z = d
            for idx, layer in enumerate(layers):
                if idx == len(layers) - 1:
                    # this is a output layer
                    prediction = layer.get_z(z)
                    prediction[prediction >= 1.0] = 1.0 - np.finfo(float).eps  # to avoid nan showing up
                    prediction[prediction <= 0.0] = 0.0 + np.finfo(float).eps
                    l1p = -l * np.log(prediction)
                    l2p = -(1.0 - l) * np.log((1.0 - prediction))
                    lcost = np.sum(l1p + l2p)
                    cost += lcost * (1.0 / float(N))

                else:

                    z = layer.get_z(z)
        if display:
            print 'cost', cost + reg
        return cost + reg

  def _testZeroDensity(self, alpha):
    """Zero isn't in the support of the gamma distribution.

    But quantized floating point math has its limits.
    TODO(bjp): Implement log-gamma sampler for small-shape distributions.

    Args:
      alpha: float shape value to test
    """
    try:
      from scipy import stats  # pylint: disable=g-import-not-at-top
    except ImportError as e:
      tf_logging.warn("Cannot test zero density proportions: %s" % e)
      return
    allowable_zeros = {
        dtypes.float16: stats.gamma(alpha).cdf(np.finfo(np.float16).tiny),
        dtypes.float32: stats.gamma(alpha).cdf(np.finfo(np.float32).tiny),
        dtypes.float64: stats.gamma(alpha).cdf(np.finfo(np.float64).tiny)
    }
    failures = []
    for use_gpu in [False, True]:
      for dt in dtypes.float16, dtypes.float32, dtypes.float64:
        sampler = self._Sampler(
            10000, alpha, 1.0, dt, use_gpu=use_gpu, seed=12345)
        x = sampler()
        allowable = allowable_zeros[dt] * x.size
        allowable = allowable * 2 if allowable < 10 else allowable * 1.05
        if np.sum(x <= 0) > allowable:
          failures += [(use_gpu, dt)]
      self.assertEqual([], failures)

    def test_exp(self):
        source = np.array([0, 0, 0, 1.0])
        result = quaternion.exp(source)
        expected = np.array([0, 0, 0, np.exp(1)])
        np.testing.assert_almost_equal(result, expected)

        source = quaternion.create_from_eulers([np.pi, 0, 0])
        result = quaternion.exp(source)
        expected = np.array([0.84147098, 0, 0, 0.54030231])
        np.testing.assert_almost_equal(result, expected)

        # Tests from the boost::math::quaternion
        source = np.array([4 * np.arctan(1), 0, 0, 0])
        result = quaternion.exp(source) + [0, 0, 0, 1.0]
        result = np.linalg.norm(result)
        expected = 2 * np.finfo(result.dtype).eps
        np.testing.assert_almost_equal(result, expected)

        source = np.array([0, 4 * np.arctan(1), 0, 0])
        result = quaternion.exp(source) + [0, 0, 0, 1.0]
        result = np.linalg.norm(result)
        expected = 2 * np.finfo(result.dtype).eps
        np.testing.assert_almost_equal(result, expected)

        source = np.array([0, 0, 4 * np.arctan(1), 0])
        result = quaternion.exp(source) + [0, 0, 0, 1.0]
        result = np.linalg.norm(result)
        expected = 2 * np.finfo(result.dtype).eps
        np.testing.assert_almost_equal(result, expected)

def elop(X, Y, op):
    """
    Compute element-wise operation of matrix :param:`X` and matrix :param:`Y`.
    
    :param X: First input matrix.
    :type X: :class:`scipy.sparse` of format csr, csc, coo, bsr, dok, lil, dia or :class:`numpy.matrix`
    :param Y: Second input matrix.
    :type Y: :class:`scipy.sparse` of format csr, csc, coo, bsr, dok, lil, dia or :class:`numpy.matrix`
    :param op: Operation to be performed. 
    :type op: `func` 
    """
    try:
        zp1 = op(0, 1) if sp.isspmatrix(X) else op(1, 0)
        zp2 = op(0, 0) 
        zp = zp1 != 0 or zp2 != 0
    except:
        zp = 0
    if sp.isspmatrix(X) or sp.isspmatrix(Y):
        return _op_spmatrix(X, Y, op) if not zp else _op_matrix(X, Y, op)
    else:
        try:
            X[X == 0] = np.finfo(X.dtype).eps
            Y[Y == 0] = np.finfo(Y.dtype).eps
        except ValueError:
            return op(np.asmatrix(X), np.asmatrix(Y))
        return op(np.asmatrix(X), np.asmatrix(Y))