add cond

Author: nschloe

.github/workflows/ci.yml

@@ -34,6 +34,6 @@ jobs:
       - name: Test with tox
         run: |
           pip install tox
-          tox
+          tox -- --cov accupy --cov-report xml --cov-report term
       - uses: codecov/codecov-action@v1
         if: ${{ matrix.python-version == '3.9' }}

README.md

@@ -36,6 +36,13 @@ p, exact, cond = accupy.generate_ill_conditioned_sum(100, 1.0e20)
 which, given a length and a target condition number, will produce an array of
 floating point numbers that is hard to sum up.
 
+Given one or two vectors, accupy can compute the condition of the sum or dot product via
+
+```python
+accupy.cond(x)
+accupy.cond(x, y)
+```
+
 accupy has the following methods for summation:
 
 - `accupy.kahan_sum(p)`: [Kahan

accupy/__init__.py

@@ -1,8 +1,13 @@
 from .dot import fdot, kdot
-from .ill_cond import generate_ill_conditioned_dot_product, generate_ill_conditioned_sum
+from .ill_cond import (
+    cond,
+    generate_ill_conditioned_dot_product,
+    generate_ill_conditioned_sum,
+)
 from .sums import decker_sum, distill, fsum, kahan_sum, knuth_sum, ksum
 
 __all__ = [
+    "cond",
     "kdot",
     "fdot",
     "generate_ill_conditioned_sum",

accupy/ill_cond.py

@@ -1,39 +1,46 @@
 import math
-from typing import Tuple
+from typing import Optional, Tuple
 
 import numpy as np
 import pyfma
-from mpmath import mp
+from numpy.typing import ArrayLike
 
+from .dot import fdot, fsum
 
-def generate_ill_conditioned_sum(
-    n: int, c: float, dps: int = 100
-) -> Tuple[np.ndarray, float, float]:
+
+def cond(x: ArrayLike, y: Optional[ArrayLike] = None) -> float:
+    """Compute the condition number of a sum (if only x is given) or a dot-product (if
+    both x and y are given).
+    """
+    if y is None:
+        return fsum(np.abs(x)) / np.abs(fsum(x))
+
+    return 2 * fdot(np.abs(x), np.abs(y)) / abs(fdot(x, y))
+
+
+def generate_ill_conditioned_sum(n: int, c: float) -> Tuple[np.ndarray, float, float]:
     # From <https://doi.org/10.1137/030601818>:
     # Ill-conditioned sums of length 2n are generated from dot products of
     # length n using Algorithm 3.3 (TwoProduct) and randomly permuting the
     # summands.
-    x, y, _, C = generate_ill_conditioned_dot_product(n, c, dps)
+    x, y, _, C = generate_ill_conditioned_dot_product(n, c)
 
     prod = x * y
     err = pyfma.fma(x, y, -prod)
     res = np.array([prod, err])
 
     out = np.random.permutation(res.flatten())
 
-    def sum_exact(p):
-        mp.dps = dps
-        return mp.fsum(p)
+    exact = fsum(out)
 
-    exact = sum_exact(out)
-    # cond = sum_exact(np.abs(out)) / abs(exact)
-    cond = C / 2
+    # condition = fsum(np.abs(out)) / abs(exact)
+    condition = C / 2
 
-    return out, exact, cond
+    return out, exact, condition
 
 
 def generate_ill_conditioned_dot_product(
-    n: int, c: float, dps: int = 100
+    n: int, c: float
 ) -> Tuple[np.ndarray, np.ndarray, float, float]:
     """n ... length of vector
     c ... target condition number
@@ -60,12 +67,6 @@ def generate_ill_conditioned_dot_product(
     x[:n2] = (2 * rx - 1) * 2 ** e
     y[:n2] = (2 * ry - 1) * 2 ** e
 
-    def dot_exact(x, y):
-        mp.dps = dps
-        # convert to list first, see
-        # <https://github.com/fredrik-johansson/mpmath/pull/385>
-        return mp.fdot(x.tolist(), y.tolist())
-
     # for i=n2+1:n and v=1:i,
     #     generate x_i, y_i such that (*) x(v)’*y(v) ~ 2^e(i-n2)
     # generate exponents for second half
@@ -76,13 +77,13 @@ def dot_exact(x, y):
         x[i] = (2 * rx[i - n2] - 1) * 2 ** e[i - n2]
         # y_i according to (*)
         y[i] = (
-            (2 * ry[i - n2] - 1) * 2 ** e[i - n2] - dot_exact(x[: i + 1], y[: i + 1])
+            (2 * ry[i - n2] - 1) * 2 ** e[i - n2] - fdot(x[: i + 1], y[: i + 1])
         ) / x[i]
 
     x, y = np.random.permutation((x, y))
     # the true dot product rounded to nearest floating point
-    d = dot_exact(x, y)
+    d = fdot(x, y)
     # the actual condition number
-    C = 2 * dot_exact(abs(x), abs(y)) / abs(d)
+    C = 2 * fdot(abs(x), abs(y)) / abs(d)
 
     return x, y, d, C

justfile

@@ -13,7 +13,7 @@ upload: clean
 	@if [ "$(git rev-parse --abbrev-ref HEAD)" != "main" ]; then exit 1; fi
 	# https://stackoverflow.com/a/58756491/353337
 	python3 -m build --sdist --wheel .
-	twine upload dist/*
+	twine upload dist/*.tar.gz
 
 publish: tag upload
 

tests/test_dot.py

@@ -1,12 +1,19 @@
 import dufte
 import matplotlib.pyplot as plt
-import numpy
+import numpy as np
 import perfplot
 import pytest
 
 import accupy
 
 
+def test_cond():
+    cond = accupy.cond([np.pi, np.e], [23225 / 8544, -355 / 113])
+    print(cond)
+    ref = 4.852507317687677e7
+    assert abs(cond - ref) < 1.0e-13 * abs(ref)
+
+
 @pytest.mark.parametrize("cond", [1.0, 1.0e15])
 def test_kdot2(cond):
     x, y, ref, _ = accupy.generate_ill_conditioned_dot_product(100, cond)
@@ -32,22 +39,22 @@ def test_accuracy_comparison_illcond(target_cond=None):
         target_cond = [10 ** k for k in range(2)]
 
     kernels = [
-        numpy.dot,
+        np.dot,
         lambda x, y: accupy.kdot(x, y, K=2),
         lambda x, y: accupy.kdot(x, y, K=3),
         accupy.fdot,
     ]
-    labels = ["numpy.dot", "accupy.kdot[2]", "accupy.kdot[3]", "accupy.fdot"]
-    data = numpy.empty((len(target_cond), len(kernels)))
-    condition_numbers = numpy.empty(len(target_cond))
-    numpy.random.seed(0)
+    labels = ["np.dot", "accupy.kdot[2]", "accupy.kdot[3]", "accupy.fdot"]
+    data = np.empty((len(target_cond), len(kernels)))
+    condition_numbers = np.empty(len(target_cond))
+    np.random.seed(0)
     for k, tc in enumerate(target_cond):
         x, y, ref, C = accupy.generate_ill_conditioned_dot_product(1000, tc)
         condition_numbers[k] = C
         data[k] = [abs(kernel(x, y) - ref) / abs(ref) for kernel in kernels]
 
     # sort
-    s = numpy.argsort(condition_numbers)
+    s = np.argsort(condition_numbers)
     condition_numbers = condition_numbers[s]
     data = data[s]
 
@@ -65,16 +72,16 @@ def test_speed_comparison1(n_range=None):
     if n_range is None:
         n_range = [2 ** k for k in range(2)]
 
-    numpy.random.seed(0)
+    np.random.seed(0)
     perfplot.plot(
-        setup=lambda n: (numpy.random.rand(n, 100), numpy.random.rand(100, n)),
+        setup=lambda n: (np.random.rand(n, 100), np.random.rand(100, n)),
         kernels=[
-            lambda xy: numpy.dot(*xy),
+            lambda xy: np.dot(*xy),
             lambda xy: accupy.kdot(*xy, K=2),
             lambda xy: accupy.kdot(*xy, K=3),
             lambda xy: accupy.fdot(*xy),
         ],
-        labels=["numpy.dot", "accupy.kdot[2]", "accupy.kdot[3]", "accupy.fdot"],
+        labels=["np.dot", "accupy.kdot[2]", "accupy.kdot[3]", "accupy.fdot"],
         n_range=n_range,
         xlabel="n",
     )
@@ -85,16 +92,16 @@ def test_speed_comparison2(n_range=None):
     if n_range is None:
         n_range = [2 ** k for k in range(2)]
 
-    numpy.random.seed(0)
+    np.random.seed(0)
     perfplot.plot(
-        setup=lambda n: (numpy.random.rand(100, n), numpy.random.rand(n, 100)),
+        setup=lambda n: (np.random.rand(100, n), np.random.rand(n, 100)),
         kernels=[
-            lambda xy: numpy.dot(*xy),
+            lambda xy: np.dot(*xy),
             lambda xy: accupy.kdot(*xy, K=2),
             lambda xy: accupy.kdot(*xy, K=3),
             lambda xy: accupy.fdot(*xy),
         ],
-        labels=["numpy.dot", "accupy.kdot[2]", "accupy.kdot[3]", "accupy.fdot"],
+        labels=["np.dot", "accupy.kdot[2]", "accupy.kdot[3]", "accupy.fdot"],
         n_range=n_range,
         xlabel="n",
         logx=True,
@@ -104,8 +111,8 @@ def test_speed_comparison2(n_range=None):
 
 
 def test_discontiguous():
-    x = numpy.random.rand(3, 10)
-    y = numpy.random.rand(3, 10)
+    x = np.random.rand(3, 10)
+    y = np.random.rand(3, 10)
     accupy.kdot(x.T, y)
     accupy.fdot(x.T, y)
 

tests/test_sums.py

@@ -1,12 +1,18 @@
 import dufte
 import matplotlib.pyplot as plt
-import numpy
+import numpy as np
 import perfplot
 import pytest
 
 import accupy
 
 
+def test_cond():
+    cond = accupy.cond([1.0, 1.0e-16, -1.0])
+    ref = 2.0e16
+    assert abs(cond - ref) < 1.0e-13 * abs(ref)
+
+
 @pytest.mark.parametrize("cond", [1.0, 1.0e15])
 def test_ksum2(cond):
     p, ref, _ = accupy.generate_ill_conditioned_sum(100, cond)
@@ -33,32 +39,32 @@ def test_accuracy_comparison_illcond(target_conds=None):
 
     kernels = [
         sum,
-        numpy.sum,
+        np.sum,
         accupy.kahan_sum,
         lambda p: accupy.ksum(p, K=2),
         lambda p: accupy.ksum(p, K=3),
         accupy.fsum,
     ]
     labels = [
         "sum",
-        "numpy.sum",
+        "np.sum",
         "accupy.kahan_sum",
         "accupy.ksum[2]",
         "accupy.ksum[3]",
         "accupy.fsum",
     ]
     colors = plt.rcParams["axes.prop_cycle"].by_key()["color"][: len(labels)]
 
-    data = numpy.empty((len(target_conds), len(kernels)))
-    condition_numbers = numpy.empty(len(target_conds))
-    numpy.random.seed(0)
+    data = np.empty((len(target_conds), len(kernels)))
+    condition_numbers = np.empty(len(target_conds))
+    np.random.seed(0)
     for k, target_cond in enumerate(target_conds):
         p, ref, C = accupy.generate_ill_conditioned_sum(1000, target_cond)
         condition_numbers[k] = C
         data[k] = [abs(kernel(p) - ref) / abs(ref) for kernel in kernels]
 
     # sort
-    s = numpy.argsort(condition_numbers)
+    s = np.argsort(condition_numbers)
     condition_numbers = condition_numbers[s]
     data = data[s]
 
@@ -77,20 +83,20 @@ def test_speed_comparison1(n_range=None):
     if n_range is None:
         n_range = [2 ** k for k in range(2)]
 
-    numpy.random.seed(0)
+    np.random.seed(0)
     perfplot.plot(
-        setup=lambda n: numpy.random.rand(n, 100),
+        setup=lambda n: np.random.rand(n, 100),
         kernels=[
             sum,
-            lambda p: numpy.sum(p, axis=0),
+            lambda p: np.sum(p, axis=0),
             accupy.kahan_sum,
             lambda p: accupy.ksum(p, K=2),
             lambda p: accupy.ksum(p, K=3),
             accupy.fsum,
         ],
         labels=[
             "sum",
-            "numpy.sum",
+            "np.sum",
             "accupy.kahan_sum",
             "accupy.ksum[2]",
             "accupy.ksum[3]",
@@ -108,20 +114,20 @@ def test_speed_comparison2(n_range=None):
     if n_range is None:
         n_range = [2 ** k for k in range(2)]
 
-    numpy.random.seed(0)
+    np.random.seed(0)
     perfplot.plot(
-        setup=lambda n: numpy.random.rand(100, n),
+        setup=lambda n: np.random.rand(100, n),
         kernels=[
             sum,
-            lambda p: numpy.sum(p, axis=0),
+            lambda p: np.sum(p, axis=0),
             accupy.kahan_sum,
             lambda p: accupy.ksum(p, K=2),
             lambda p: accupy.ksum(p, K=3),
             accupy.fsum,
         ],
         labels=[
             "sum",
-            "numpy.sum",
+            "np.sum",
             "accupy.kahan_sum",
             "accupy.ksum[2]",
             "accupy.ksum[3]",
@@ -134,25 +140,25 @@ def test_speed_comparison2(n_range=None):
 
 
 def test_knuth_sum():
-    a16 = numpy.float16(1.0e1)
-    b16 = numpy.float16(1.0e-1)
+    a16 = np.float16(1.0e1)
+    b16 = np.float16(1.0e-1)
     x16, y16 = accupy.knuth_sum(a16, b16)
-    xy = numpy.float64(x16) + numpy.float64(y16)
-    ab = numpy.float64(a16) + numpy.float64(b16)
+    xy = np.float64(x16) + np.float64(y16)
+    ab = np.float64(a16) + np.float64(b16)
     assert abs(xy - ab) < 1.0e-15 * ab
 
 
 def test_decker_sum():
-    a16 = numpy.float16(1.0e1)
-    b16 = numpy.float16(1.0e-1)
+    a16 = np.float16(1.0e1)
+    b16 = np.float16(1.0e-1)
     x16, y16 = accupy.decker_sum(a16, b16)
-    xy = numpy.float64(x16) + numpy.float64(y16)
-    ab = numpy.float64(a16) + numpy.float64(b16)
+    xy = np.float64(x16) + np.float64(y16)
+    ab = np.float64(a16) + np.float64(b16)
     assert abs(xy - ab) < 1.0e-15 * ab
 
 
 def test_discontiguous():
-    x = numpy.random.rand(3, 10).T
+    x = np.random.rand(3, 10).T
     accupy.ksum(x.T)
     accupy.fsum(x.T)
 

tox.ini

@@ -9,4 +9,4 @@ deps =
     pytest
     pytest-cov
 commands =
-    pytest --cov {envsitepackagesdir}/accupy --cov-report xml --cov-report term
+    pytest {posargs}