Python numpy.nextafter函数代码示例

def make_strictly_feasible(x, lb, ub, rstep=1e-10):
    """Shift a point to the interior of a feasible region.
    
    Each element of the returned vector is at least at a relative distance
    `rstep` from the closest bound. If ``rstep=0`` then `np.nextafter` is used.
    """
    x_new = x.copy()

    active = find_active_constraints(x, lb, ub, rstep)
    lower_mask = np.equal(active, -1)
    upper_mask = np.equal(active, 1)

    if rstep == 0:
        x_new[lower_mask] = np.nextafter(lb[lower_mask], ub[lower_mask])
        x_new[upper_mask] = np.nextafter(ub[upper_mask], lb[upper_mask])
    else:
        x_new[lower_mask] = (lb[lower_mask] +
                             rstep * np.maximum(1, np.abs(lb[lower_mask])))
        x_new[upper_mask] = (ub[upper_mask] -
                             rstep * np.maximum(1, np.abs(ub[upper_mask])))

    tight_bounds = (x_new < lb) | (x_new > ub)
    x_new[tight_bounds] = 0.5 * (lb[tight_bounds] + ub[tight_bounds])

    return x_new

def make_strictly_feasible(x, lb, ub, rstep=0):
    """Shift the point in the slightest possible way to the interior.

    If ``rstep=0`` the function uses np.nextafter, otherwise `rstep` is
    multiplied by absolute value of the bound.

    The utility of this function is questionable to me. Maybe bigger shifts
    should be used, or maybe this function is not necessary at all despite
    theoretical requirement of our interior point algorithm.
    """
    x_new = x.copy()

    m = x <= lb
    if rstep == 0:
        x_new[m] = np.nextafter(lb[m], ub[m])
    else:
        x_new[m] = lb[m] + rstep * (1 + np.abs(lb[m]))

    m = x >= ub
    if rstep == 0:
        x_new[m] = np.nextafter(ub[m], lb[m])
    else:
        x_new[m] = ub[m] - rstep * (1 + np.abs(ub[m]))

    return x_new

def nextafter(x, direction, dtype, itemsize):
    """Return the next representable neighbor of x in the appropriate
    direction."""

    assert direction in [-1, 0, +1]
    assert dtype.kind == "S" or type(x) in (bool, int, int, float)

    if direction == 0:
        return x

    if dtype.kind == "S":
        return string_next_after(x, direction, itemsize)

    if dtype.kind in ['b']:
        return bool_type_next_after(x, direction, itemsize)
    elif dtype.kind in ['i', 'u']:
        return int_type_next_after(x, direction, itemsize)
    elif dtype.kind == "f":
        if direction < 0:
            return numpy.nextafter(x, x - 1)
        else:
            return numpy.nextafter(x, x + 1)

    # elif dtype.name == "float32":
    #    if direction < 0:
    #        return PyNextAfterF(x,x-1)
    #    else:
    #        return PyNextAfterF(x,x + 1)
    # elif dtype.name == "float64":
    #    if direction < 0:
    #        return PyNextAfter(x,x-1)
    #    else:
    #        return PyNextAfter(x,x + 1)

    raise TypeError("data type ``%s`` is not supported" % dtype)

    def __init__(self, value_1, value_2=None):
        # use Decimal as exact value holder (because of arbitrary precision)
        from decimal import Decimal
        # nextafter(x, y) returns next machine number after x in direction of y
        from numpy import nextafter

        # creating interval from middle value
        if not value_2:
            exact = Decimal(value_1)
            float_repr = Decimal("{0:0.70f}".format(float(exact)))
            if exact == float_repr:
                self.lv = float(float_repr)
                self.rv = float(float_repr)
            elif exact > float_repr:
                self.lv = float(float_repr)
                self.rv = nextafter(self.lv, float('Inf'))
            elif exact < float_repr:
                self.rv = float(float_repr)
                self.lv = nextafter(self.rv, -float('Inf'))
        # creating interval from left and right edge
        else:
            exact_left = Decimal(value_1)
            exact_right = Decimal(value_2)
            if exact_left > exact_right:
                exact_left, exact_right = exact_right, exact_left
            float_repr_left = Decimal(float(exact_left))
            float_repr_right = Decimal(float(exact_right))
            if exact_left < float_repr_left:
                self.lv = nextafter(float(float_repr_left), -float('Inf'))
            else:
                self.lv = float(float_repr_left)
            if exact_right > float_repr_right:
                self.rv = nextafter(float(float_repr_right), float('Inf'))
            else:
                self.rv = float(float_repr_right)

 def _logpmf(self, x, mu, alpha, p):
     mu_p = mu ** (p - 1.)
     a1 = np.maximum(np.nextafter(0, 1), 1 + alpha * mu_p)
     a2 = np.maximum(np.nextafter(0, 1), mu + (a1 - 1.) * x)
     logpmf_ = np.log(mu) + (x - 1.) * np.log(a2)
     logpmf_ -=  x * np.log(a1) + gammaln(x + 1.) + a2 / a1
     return logpmf_

    def test_half_fpe(self):
        oldsettings = np.seterr(all="raise")
        try:
            sx16 = np.array((1e-4,), dtype=float16)
            bx16 = np.array((1e4,), dtype=float16)
            sy16 = float16(1e-4)
            by16 = float16(1e4)

            # Underflow errors
            assert_raises_fpe("underflow", lambda a, b: a * b, sx16, sx16)
            assert_raises_fpe("underflow", lambda a, b: a * b, sx16, sy16)
            assert_raises_fpe("underflow", lambda a, b: a * b, sy16, sx16)
            assert_raises_fpe("underflow", lambda a, b: a * b, sy16, sy16)
            assert_raises_fpe("underflow", lambda a, b: a / b, sx16, bx16)
            assert_raises_fpe("underflow", lambda a, b: a / b, sx16, by16)
            assert_raises_fpe("underflow", lambda a, b: a / b, sy16, bx16)
            assert_raises_fpe("underflow", lambda a, b: a / b, sy16, by16)
            assert_raises_fpe("underflow", lambda a, b: a / b, float16(2.0 ** -14), float16(2 ** 11))
            assert_raises_fpe("underflow", lambda a, b: a / b, float16(-2.0 ** -14), float16(2 ** 11))
            assert_raises_fpe("underflow", lambda a, b: a / b, float16(2.0 ** -14 + 2 ** -24), float16(2))
            assert_raises_fpe("underflow", lambda a, b: a / b, float16(-2.0 ** -14 - 2 ** -24), float16(2))
            assert_raises_fpe("underflow", lambda a, b: a / b, float16(2.0 ** -14 + 2 ** -23), float16(4))

            # Overflow errors
            assert_raises_fpe("overflow", lambda a, b: a * b, bx16, bx16)
            assert_raises_fpe("overflow", lambda a, b: a * b, bx16, by16)
            assert_raises_fpe("overflow", lambda a, b: a * b, by16, bx16)
            assert_raises_fpe("overflow", lambda a, b: a * b, by16, by16)
            assert_raises_fpe("overflow", lambda a, b: a / b, bx16, sx16)
            assert_raises_fpe("overflow", lambda a, b: a / b, bx16, sy16)
            assert_raises_fpe("overflow", lambda a, b: a / b, by16, sx16)
            assert_raises_fpe("overflow", lambda a, b: a / b, by16, sy16)
            assert_raises_fpe("overflow", lambda a, b: a + b, float16(65504), float16(17))
            assert_raises_fpe("overflow", lambda a, b: a - b, float16(-65504), float16(17))
            assert_raises_fpe("overflow", np.nextafter, float16(65504), float16(np.inf))
            assert_raises_fpe("overflow", np.nextafter, float16(-65504), float16(-np.inf))
            assert_raises_fpe("overflow", np.spacing, float16(65504))

            # Invalid value errors
            assert_raises_fpe("invalid", np.divide, float16(np.inf), float16(np.inf))
            assert_raises_fpe("invalid", np.spacing, float16(np.inf))
            assert_raises_fpe("invalid", np.spacing, float16(np.nan))
            assert_raises_fpe("invalid", np.nextafter, float16(np.inf), float16(0))
            assert_raises_fpe("invalid", np.nextafter, float16(-np.inf), float16(0))
            assert_raises_fpe("invalid", np.nextafter, float16(0), float16(np.nan))

            # These should not raise
            float16(65472) + float16(32)
            float16(2 ** -13) / float16(2)
            float16(2 ** -14) / float16(2 ** 10)
            np.spacing(float16(-65504))
            np.nextafter(float16(65504), float16(-np.inf))
            np.nextafter(float16(-65504), float16(np.inf))
            float16(2 ** -14) / float16(2 ** 10)
            float16(-2 ** -14) / float16(2 ** 10)
            float16(2 ** -14 + 2 ** -23) / float16(2)
            float16(-2 ** -14 - 2 ** -23) / float16(2)
        finally:
            np.seterr(**oldsettings)

    def forward_cpu(self, inputs):
        U, points = inputs
        batch_size, height, width = U.shape

        # Points just on the boundary are slightly (i.e. nextafter in float32)
        # moved inward to simplify the implementation
        points = points.copy()
        on_boundary = (points == 0)
        points[on_boundary] = np.nextafter(points[on_boundary], np.float32(1))
        x = points[:, 0]
        y = points[:, 1]
        on_boundary = (x == (width - 1))
        x[on_boundary] = np.nextafter(x[on_boundary], np.float32(0))
        on_boundary = (y == (height - 1))
        y[on_boundary] = np.nextafter(y[on_boundary], np.float32(0))

        batch_axis = np.expand_dims(np.arange(batch_size), 1)
        points_floor = np.floor(points)
        x_l = points_floor[:, 0].astype(np.int32)
        y_l = points_floor[:, 1].astype(np.int32)
        x_l = np.clip(x_l, 0, width - 1)
        y_l = np.clip(y_l, 0, height - 1)
        x_h = np.clip(x_l + 1, 0, width - 1)
        y_h = np.clip(y_l + 1, 0, height - 1)

        weight = 1.0 - (points - points_floor)
        weight_x_l = weight[:, 0]
        weight_y_l = weight[:, 1]
        weight_x_h = 1 - weight_x_l
        weight_y_h = 1 - weight_y_l

        # remove points outside of the (source) image region
        # by setting their weights to 0
        x_invalid = np.logical_or(x < 0, (width - 1) < x)
        y_invalid = np.logical_or(y < 0, (height - 1) < y)
        invalid = np.logical_or(x_invalid, y_invalid)
        weight_x_l[invalid] = 0
        weight_y_l[invalid] = 0
        weight_x_h[invalid] = 0
        weight_y_h[invalid] = 0

        U_y_l = (weight_x_l * U[batch_axis, y_l, x_l] +
                 weight_x_h * U[batch_axis, y_l, x_h])
        U_y_h = (weight_x_l * U[batch_axis, y_h, x_l] +
                 weight_x_h * U[batch_axis, y_h, x_h])
        V = weight_y_l * U_y_l + weight_y_h * U_y_h

        self.x_l = x_l
        self.y_l = y_l
        self.x_h = x_h
        self.y_h = y_h
        self.weight_x_l = weight_x_l
        self.weight_y_l = weight_y_l
        self.weight_x_h = weight_x_h
        self.weight_y_h = weight_y_h
        return (V,)

def _test_nextafter(t):
    one = t(1)
    two = t(2)
    zero = t(0)
    eps = np.finfo(t).eps
    assert_(np.nextafter(one, two) - one == eps)
    assert_(np.nextafter(one, zero) - one < 0)
    assert_(np.isnan(np.nextafter(np.nan, one)))
    assert_(np.isnan(np.nextafter(one, np.nan)))
    assert_(np.nextafter(one, one) == one)

def test_nextafter():
    for t in [np.float32, np.float64, np.longdouble]:
        one = t(1)
        two = t(2)
        zero = t(0)
        eps = np.finfo(t).eps
        assert np.nextafter(one, two) - one == eps
        assert np.nextafter(one, zero) - one < 0
        assert np.isnan(np.nextafter(np.nan, one))
        assert np.isnan(np.nextafter(one, np.nan))
        assert np.nextafter(one, one) == one

def ranges_to_weight_table(ranges):
    """
    Create a table of weights from ranges. Include only edge points of ranges.
    Include each edge point twice: once as values within the range and zero
    value outside the range (with this output the weights can easily be interpolated).

    Weights of overlapping intervals are summed.

    Assumes 64-bit floats.

    :param ranges: list of triples (edge1, edge2, weight)
    :return: an Orange.data.Table
    """

    values = {}

    inf = float("inf")
    minf = float("-inf")

    def dict_to_numpy(d):
        x = []
        y = []
        for a, b in d.items():
            x.append(a)
            y.append(b)
        return np.array(x), np.array([y])

    for l, r, w in ranges:
        l, r = min(l, r), max(l, r)
        positions = [nextafter(l, minf), l, r, nextafter(r, inf)]
        weights = [0., float(w), float(w), 0.]

        all_positions = list(set(positions) | set(values))  # new and old positions

        # current values on all position
        x, y = dict_to_numpy(values)
        current = interp1d_with_unknowns_numpy(x, y, all_positions)[0]
        current[np.isnan(current)] = 0

        # new values on all positions
        new = interp1d_with_unknowns_numpy(np.array(positions), np.array([weights]),
                                           all_positions)[0]
        new[np.isnan(new)] = 0

        # update values
        for p, f in zip(all_positions, current + new):
            values[p] = f

    x, y = dict_to_numpy(values)
    dom = Orange.data.Domain([Orange.data.ContinuousVariable(name=str(float(a))) for a in x])
    data = Orange.data.Table.from_numpy(dom, y)
    return data

    def test_half_conversion_rounding(self, float_t, shift, offset):
        # Assumes that round to even is used during casting.
        max_pattern = np.float16(np.finfo(np.float16).max).view(np.uint16)

        # Test all (positive) finite numbers, denormals are most interesting
        # however:
        f16s_patterns = np.arange(0, max_pattern+1, dtype=np.uint16)
        f16s_float = f16s_patterns.view(np.float16).astype(float_t)

        # Shift the values by half a bit up or a down (or do not shift),
        if shift == "up":
            f16s_float = 0.5 * (f16s_float[:-1] + f16s_float[1:])[1:]
        elif shift == "down":
            f16s_float = 0.5 * (f16s_float[:-1] + f16s_float[1:])[:-1]
        else:
            f16s_float = f16s_float[1:-1]

        # Increase the float by a minimal value:
        if offset == "up":
            f16s_float = np.nextafter(f16s_float, float_t(1e50))
        elif offset == "down":
            f16s_float = np.nextafter(f16s_float, float_t(-1e50))

        # Convert back to float16 and its bit pattern:
        res_patterns = f16s_float.astype(np.float16).view(np.uint16)

        # The above calculations tries the original values, or the exact
        # mid points between the float16 values. It then further offsets them
        # by as little as possible. If no offset occurs, "round to even"
        # logic will be necessary, an arbitrarily small offset should cause
        # normal up/down rounding always.

        # Calculate the expecte pattern:
        cmp_patterns = f16s_patterns[1:-1].copy()

        if shift == "down" and offset != "up":
            shift_pattern = -1
        elif shift == "up" and offset != "down":
            shift_pattern = 1
        else:
            # There cannot be a shift, either shift is None, so all rounding
            # will go back to original, or shift is reduced by offset too much.
            shift_pattern = 0

        # If rounding occurs, is it normal rounding or round to even?
        if offset is None:
            # Round to even occurs, modify only non-even, cast to allow + (-1)
            cmp_patterns[0::2].view(np.int16)[...] += shift_pattern
        else:
            cmp_patterns.view(np.int16)[...] += shift_pattern

        assert_equal(res_patterns, cmp_patterns)

 def compute_likelihoods(self, PLCs, FLCs):
     K = self.K()
     N = self.N()
     future_given_state_probs = np.nextafter(self.f_hat_conditional_densities(FLCs, label="PDF_FLCs"), 1.)
     state_given_past_probs = np.nextafter(np.vstack([self.PLC_densities(j, PLCs) for j in range(K)]), 1.).T        
     ''' Weight by state likelihood '''
     n_hats = self.W.sum(axis=0) / N
     state_given_past_probs *= n_hats
     state_given_past_probs = np.nextafter(state_given_past_probs, 1.)
     ''' Normalize '''
     state_given_past_probs /= np.expand_dims(np.sum(state_given_past_probs, axis=1), axis=1)
     ''' Return mixed likelihoods '''
     return np.nextafter(np.sum(np.multiply(state_given_past_probs, future_given_state_probs), axis=1), 1.)

    def testBearingToValueOnEquator(self):
        """Test if bearingTo() returns the expected value from a point on the equator
        """
        lon0 = 90.0
        lat0 = 0.0   # These tests only work from the equator.
        arcLen = 10.0

        trials = [
            # Along celestial equator
            dict(lon=lon0, lat=lat0, bearing=0.0,
                 lonEnd=lon0+arcLen, latEnd=lat0),
            # Along a meridian
            dict(lon=lon0, lat=lat0, bearing=90.0,
                 lonEnd=lon0, latEnd=lat0+arcLen),
            # 180 degree arc (should go to antipodal point)
            dict(lon=lon0, lat=lat0, bearing=45.0,
                 lonEnd=lon0+180.0, latEnd=-lat0),
            #
            dict(lon=lon0, lat=lat0, bearing=45.0,
                 lonEnd=lon0+90.0, latEnd=lat0 + 45.0),
            dict(lon=lon0, lat=lat0, bearing=225.0,
                 lonEnd=lon0-90.0, latEnd=lat0 - 45.0),
            dict(lon=lon0, lat=np.nextafter(-90.0, inf),
                 bearing=90.0, lonEnd=lon0, latEnd=0.0),
            dict(lon=lon0, lat=np.nextafter(-90.0, inf),
                 bearing=0.0, lonEnd=lon0 + 90.0, latEnd=0.0),
            # Argument at a pole should work
            dict(lon=lon0, lat=lat0, bearing=270.0, lonEnd=lon0, latEnd=-90.0),
            # Support for non-finite values
            dict(lon=lon0, lat=nan, bearing=nan, lonEnd=lon0, latEnd=45.0),
            dict(lon=lon0, lat=lat0, bearing=nan, lonEnd=nan, latEnd=90.0),
            dict(lon=inf, lat=lat0, bearing=nan, lonEnd=lon0, latEnd=42.0),
            dict(lon=lon0, lat=lat0, bearing=nan, lonEnd=-inf, latEnd=42.0),
        ]

        for trial in trials:
            origin = SpherePoint(trial['lon']*degrees, trial['lat']*degrees)
            end = SpherePoint(trial['lonEnd']*degrees, trial['latEnd']*degrees)
            bearing = origin.bearingTo(end)

            self.assertIsInstance(bearing, geom.Angle)
            if origin.isFinite() and end.isFinite():
                self.assertGreaterEqual(bearing.asDegrees(), 0.0)
                self.assertLess(bearing.asDegrees(), 360.0)
            if origin.separation(end).asDegrees() != 180.0:
                if not math.isnan(trial['bearing']):
                    self.assertAlmostEqual(
                        trial['bearing'], bearing.asDegrees(), 12)
                else:
                    self.assertTrue(math.isnan(bearing.asRadians()))

def similarity(v1, v2):
    # v1 and v2 are vectors
    eps = np.nextafter(0, 1)  # smallest float above zero
    dot = np.dot(v1, v2)
    dot /= max(npext.norm(v1), eps)
    dot /= max(npext.norm(v2), eps)
    return dot

def _compute_lwork(routine, *args, **kwargs):
    """
    Round floating-point lwork returned by lapack to integer.

    Several LAPACK routines compute optimal values for LWORK, which
    they return in a floating-point variable. However, for large
    values of LWORK, single-precision floating point is not sufficient
    to hold the exact value --- some LAPACK versions (<= 3.5.0 at
    least) truncate the returned integer to single precision and in
    some cases this can be smaller than the required value.
    """
    wi = routine(*args, **kwargs)
    if len(wi) < 2:
        raise ValueError("")
    info = wi[-1]
    if info != 0:
        raise ValueError("Internal work array size computation failed: " "%d" % (info,))

    lwork = [w.real for w in wi[:-1]]

    dtype = getattr(routine, "dtype", None)
    if dtype == _np.float32 or dtype == _np.complex64:
        # Single-precision routine -- take next fp value to work
        # around possible truncation in LAPACK code
        lwork = _np.nextafter(lwork, _np.inf, dtype=_np.float32)

    lwork = _np.array(lwork, _np.int64)
    if _np.any(_np.logical_or(lwork < 0, lwork > _np.iinfo(_np.int32).max)):
        raise ValueError(
            "Too large work array required -- computation cannot " "be performed with standard 32-bit LAPACK."
        )
    lwork = lwork.astype(_np.int32)
    if lwork.size == 1:
        return lwork[0]
    return lwork