Merge pull request #530 from qutech/issues/529_from_float

Issues/529 from float

ReleaseNotes.txt

@@ -1,3 +1,8 @@
+## pending/current ##
+
+ - General:
+    - Unify `TimeType.from_float` between fractions and gmpy2 backend behaviour.
+
 ## 0.5 ##
 
  - General:

qupulse/utils/numeric.py

@@ -0,0 +1,100 @@
+from typing import Tuple, Type
+from numbers import Rational
+from math import gcd
+
+
+def lcm(a: int, b: int):
+    """least common multiple"""
+    return a * b // gcd(a, b)
+
+
+def _approximate_int(alpha_num: int, d_num: int, den: int) -> Tuple[int, int]:
+    """Find the best fraction approximation of alpha_num / den with an error smaller d_num / den. Best means the
+    fraction with the smallest denominator.
+
+    Algorithm from https://link.springer.com/content/pdf/10.1007%2F978-3-540-72914-3.pdf
+
+    Args:s
+        alpha_num: Numerator of number to approximate. 0 < alpha_num < den
+        d_num: Numerator of allowed absolute error.
+        den: Denominator of both numbers above.
+
+    Returns:
+        (numerator, denominator)
+    """
+    assert 0 < alpha_num < den
+
+    lower_num = alpha_num - d_num
+    upper_num = alpha_num + d_num
+
+    p_a, q_a = 0, 1
+    p_b, q_b = 1, 1
+
+    p_full, q_full = p_b, q_b
+
+    to_left = True
+
+    while True:
+
+        # compute the number of steps to the left
+        x_num = den * p_b - alpha_num * q_b
+        x_den = -den * p_a + alpha_num * q_a
+        x = (x_num + x_den - 1) // x_den  # ceiling division
+
+        p_full += x * p_a
+        q_full += x * q_a
+
+        p_prev = p_full - p_a
+        q_prev = q_full - q_a
+
+        # check whether we have a valid approximation
+        if (q_full * lower_num < p_full * den < q_full * upper_num or
+                q_prev * lower_num < p_prev * den < q_prev * upper_num):
+            bound_num = upper_num if to_left else lower_num
+
+            k_num = den * p_b - bound_num * q_b
+            k_den = bound_num * q_a - den * p_a
+            k = (k_num // k_den) + 1
+
+            return p_b + k * p_a, q_b + k * q_a
+
+        # update the interval
+        p_a = p_prev
+        q_a = q_prev
+
+        p_b = p_full
+        q_b = q_full
+
+        to_left = not to_left
+
+
+def approximate_rational(x: Rational, abs_err: Rational, fraction_type: Type[Rational]) -> Rational:
+    """Return the fraction with the smallest denominator in (x - abs_err, x + abs_err)"""
+    if abs_err <= 0:
+        raise ValueError('abs_err must be > 0')
+
+    xp, xq = x.numerator, x.denominator
+    if xq == 1:
+        return x
+
+    dp, dq = abs_err.numerator, abs_err.denominator
+
+    # separate integer part. alpha_num is guaranteed to be < xq
+    n, alpha_num = divmod(xp, xq)
+
+    # find common denominator of alpha_num / xq and dp / dq
+    den = lcm(xq, dq)
+    alpha_num = alpha_num * den // xq
+    d_num = dp * den // dq
+
+    if alpha_num < d_num:
+        p, q = 0, 1
+    else:
+        p, q = _approximate_int(alpha_num, d_num, den)
+
+    return fraction_type(p + n * q, q)
+
+
+def approximate_double(x: float, abs_err: float, fraction_type: Type[Rational]) -> Rational:
+    """Return the fraction with the smallest denominator in (x - abs_err, x + abs_err)."""
+    return approximate_rational(fraction_type(x), fraction_type(abs_err), fraction_type=fraction_type)

qupulse/utils/types.py

@@ -10,6 +10,8 @@
 
 import numpy
 
+import qupulse.utils.numeric as qupulse_numeric
+
 __all__ = ["MeasurementWindow", "ChannelID", "HashableNumpyArray", "TimeType", "time_from_float", "DocStringABCMeta",
            "SingletonABCMeta", "SequenceProxy"]
 
@@ -18,6 +20,8 @@
 
 try:
     import gmpy2
+    qupulse_numeric.FractionType = gmpy2.mpq
+
 except ImportError:
     gmpy2 = None
 
@@ -182,31 +186,31 @@ def from_float(cls, value: float, absolute_error: typing.Optional[float] = None)
             absolute_error:
                 - :obj:`None`: Use `str(value)` as a proxy to get consistent precision
                 - 0: Return the exact value of the float i.e. float(0.8) == 3602879701896397 / 4503599627370496
-                - 0 < `absolute_error` <= 1: Use `absolute_error` to limit the denominator
+                - 0 < `absolute_error` <= 1: Return the best approximation to `value` within `(value - absolute_error,
+                value + absolute_error)`. The best approximation is defined as the fraction with the smallest
+                denominator.
 
         Raises:
             ValueError: If `absolute_error` is not None and not 0 <= `absolute_error` <=  1
         """
         # gmpy2 is at least an order of magnitude faster than fractions.Fraction
         if absolute_error is None:
-            # this method utilizes the 'print as many digits as necessary to destinguish between all floats'
+            # this method utilizes the 'print as many digits as necessary to distinguish between all floats'
             # functionality of str
             if type(value) in (cls, cls._InternalType, fractions.Fraction):
                 return cls(value)
             else:
-                return cls(cls._to_internal(str(value).replace('e', 'E')))
+                # .upper() is a bit faster than replace('e', 'E') which gmpy2.mpq needs
+                return cls(cls._to_internal(str(value).upper()))
 
         elif absolute_error == 0:
             return cls(cls._to_internal(value))
         elif absolute_error < 0:
-            raise ValueError('absolute_error needs to be at least 0')
+            raise ValueError('absolute_error needs to be > 0')
         elif absolute_error > 1:
-            raise ValueError('absolute_error needs to be smaller 1')
+            raise ValueError('absolute_error needs to be <= 1')
         else:
-            if cls._InternalType is fractions.Fraction:
-                return fractions.Fraction(value).limit_denominator(int(1 / absolute_error))
-            else:
-                return cls(gmpy2.f2q(value, absolute_error))
+            return cls(qupulse_numeric.approximate_double(value, absolute_error, fraction_type=cls._InternalType))
 
     @classmethod
     def from_fraction(cls, numerator: int, denominator: int) -> 'TimeType':

tests/utils/numeric_tests.py

@@ -0,0 +1,111 @@
+import unittest
+from typing import Callable, List, Iterator, Tuple
+import random
+from fractions import Fraction
+from collections import deque
+from itertools import islice
+
+from qupulse.utils.numeric import approximate_rational, approximate_double
+
+
+def stern_brocot_sequence() -> Iterator[int]:
+    sb = deque([1, 1])
+    while True:
+        sb += [sb[0] + sb[1], sb[1]]
+        yield sb.popleft()
+
+
+def stern_brocot_tree(depth: int) -> List[Fraction]:
+    """see wikipedia article"""
+    fractions = [Fraction(0), Fraction(1)]
+
+    seq = stern_brocot_sequence()
+    next(seq)
+
+    for n in range(depth):
+        for _ in range(n + 1):
+            p = next(seq)
+            q = next(seq)
+            fractions.append(Fraction(p, q))
+
+    return sorted(fractions)
+
+
+def window(iterable, n):
+    assert n > 0
+    it = iter(iterable)
+    state = deque(islice(it, 0, n - 1), maxlen=n)
+    for new_element in it:
+        state.append(new_element)
+        yield tuple(state)
+        state.popleft()
+
+
+def uniform_without_bounds(rng, a, b):
+    result = a
+    while not a < result < b:
+        result = rng.uniform(a, b)
+    return result
+
+
+def generate_test_pairs(depth: int, seed) -> Iterator[Tuple[Tuple[Fraction, Fraction], Fraction]]:
+    rng = random.Random(seed)
+    tree = stern_brocot_tree(depth)
+    extended_tree = [float('-inf')] + tree + [float('inf')]
+
+    # values map to themselves
+    for a, b, c in window(tree, 3):
+        err = min(b - a, c - b)
+        yield (b, err), b
+
+    for prev, a, b, upcom in zip(extended_tree, extended_tree[1:], extended_tree[2:], extended_tree[3:]):
+        mid = (a + b) / 2
+
+        low = Fraction(uniform_without_bounds(rng, a, mid))
+        err = min(mid - a, low - prev)
+        yield (low, err), a
+
+        high = Fraction(uniform_without_bounds(rng, mid, b))
+        err = min(b - mid, upcom - high)
+        yield (high, err), b
+
+
+class ApproximationTests(unittest.TestCase):
+    def test_approximate_rational(self):
+        """Use Stern-Brocot tree and rng to generate test cases where we know the result"""
+        depth = 70  # equivalent to 7457 test cases
+        test_pairs = list(generate_test_pairs(depth, seed=42))
+
+        for offset in (-2, -1, 0, 1, 2):
+            for (x, abs_err), result in test_pairs:
+                expected = result + offset
+                result = approximate_rational(x + offset, abs_err, Fraction)
+                self.assertEqual(expected, result)
+
+        with self.assertRaises(ValueError):
+            approximate_rational(Fraction(3, 1), Fraction(0, 100), Fraction)
+
+        with self.assertRaises(ValueError):
+            approximate_rational(Fraction(3, 1), Fraction(-1, 100), Fraction)
+
+        x = Fraction(3, 1)
+        abs_err = Fraction(1, 100)
+        self.assertIs(x, approximate_rational(x, abs_err, Fraction))
+
+    def test_approximate_double(self):
+        test_values = [
+            ((.1, .05), Fraction(1, 7)),
+            ((.12, .005), Fraction(2, 17)),
+            ((.15, .005), Fraction(2, 13)),
+            # .111_111_12, 0.000_000_005
+            ((.11111112, 0.000000005), Fraction(888890, 8000009)),
+            ((.111125, 0.0000005), Fraction(859, 7730)),
+            ((2.50000000008, .1), Fraction(5, 2))
+        ]
+
+        for (x, err), expected in test_values:
+            result = approximate_double(x, err, Fraction)
+            self.assertEqual(expected, result, msg='{x} ± {err} results in {result} '
+                                                   'which is not the expected {expected}'.format(x=x, err=err,
+                                                                                                 result=result,
+                                                                                                 expected=expected))

tests/utils/time_type_tests.py

@@ -87,6 +87,8 @@ def test_fraction_time_from_float_exact_fallback(self):
     def assert_fraction_time_from_float_with_precision_works(self, time_type):
         self.assertEqual(time_type.from_float(1000000 / 1000001, 1e-5),
                          fractions.Fraction(1))
+        self.assertEqual(time_type.from_float(2.50000000000008, absolute_error=1e-10),
+                         time_type.from_fraction(5, 2))
 
     def test_fraction_time_from_float_with_precision(self):
         self.assert_fraction_time_from_float_with_precision_works(qutypes.TimeType)
@@ -126,10 +128,10 @@ def test_from_float_no_extra_args_fallback(self):
         self.assert_from_float_exact_works(self.fallback_qutypes.TimeType)
 
     def test_from_float_exceptions(self):
-        with self.assertRaisesRegex(ValueError, 'at least 0'):
+        with self.assertRaisesRegex(ValueError, '> 0'):
             qutypes.time_from_float(.8, -1)
 
-        with self.assertRaisesRegex(ValueError, 'smaller 1'):
+        with self.assertRaisesRegex(ValueError, '<= 1'):
             qutypes.time_from_float(.8, 2)