sigma-py / accupy Goto Github PK

View Code? Open in Web Editor NEW

102.0 3.0 15.0 1.09 MB

:dart: Accurate sums and dot products for Python.

License: GNU General Public License v3.0

Python 77.92% C++ 22.08%

mathematics numerical-methods engineering accuracy numpy python pypi

accupy's Introduction

Accurate sums and (dot) products for Python.

Sums

Summing up values in a list can get tricky if the values are floating point numbers; digit cancellation can occur and the result may come out wrong. A classical example is the sum

1.0e16 + 1.0 - 1.0e16

The actual result is 1.0, but in double precision, this will result in 0.0. While in this example the failure is quite obvious, it can get a lot more tricky than that. accupy provides

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

accupy.cond(x)
accupy.cond(x, y)

accupy has the following methods for summation:

accupy.kahan_sum(p): Kahan summation
accupy.fsum(p): A vectorization wrapper around math.fsum (which uses Shewchuck's algorithm [1] (see also here)).
accupy.ksum(p, K=2): Summation in K-fold precision (from [2])

All summation methods sum the first dimension of a multidimensional NumPy array.

Let's compare them.

Accuracy comparison (sum)

As expected, the naive sum performs very badly with ill-conditioned sums; likewise for numpy.sum which uses pairwise summation. Kahan summation not significantly better; this, too, is expected.

Computing the sum with 2-fold accuracy in accupy.ksum gives the correct result if the condition is at most in the range of machine precision; further increasing K helps with worse conditions.

Shewchuck's algorithm in math.fsum always gives the correct result to full floating point precision.

Runtime comparison (sum)

We compare more and more sums of fixed size (above) and larger and larger sums, but a fixed number of them (below). In both cases, the least accurate method is the fastest (numpy.sum), and the most accurate the slowest (accupy.fsum).

Dot products

accupy has the following methods for dot products:

accupy.fdot(p): A transformation of the dot product of length n into a sum of length 2n, computed with math.fsum
accupy.kdot(p, K=2): Dot product in K-fold precision (from [2])

Let's compare them.

Accuracy comparison (dot)

accupy can construct ill-conditioned dot products with

x, y, exact, cond = accupy.generate_ill_conditioned_dot_product(100, 1.0e20)

With this, the accuracy of the different methods is compared.

As for sums, numpy.dot is the least accurate, followed by instanced of kdot. fdot is provably accurate up into the last digit

Runtime comparison (dot)

NumPy's numpy.dot is much faster than all alternatives provided by accupy. This is because the bookkeeping of truncation errors takes more steps, but mostly because of NumPy's highly optimized dot implementation.

References

Dependencies

accupy needs the C++ Eigen library, provided in Debian/Ubuntu by libeigen3-dev.

Installation

accupy is available from the Python Package Index, so with

pip install accupy

you can install.

Testing

To run the tests, just check out this repository and type

MPLBACKEND=Agg pytest

accupy's People

Contributors

Stargazers

Watchers

Forkers

kostyfisik hhc2tech toddrme2178 risicle ahmed-f-alrefaie gugar20 bjodah numsim1415 diego-hayashi tumblerxi txetdesrever

accupy's Issues

check Twofold fast summation by Evgeny Latkin

https://arxiv.org/pdf/1401.0248.pdf

test failure on macos: test_kdot2

On macos 10.14, having built accupy 0.3.2 using clang 7.1.0, I get a precision error in the test test/test_dot.py::test_kdot2[1000000000000000.0]. Output:

============================= test session starts ==============================
platform darwin -- Python 3.7.9, pytest-5.4.3, py-1.9.0, pluggy-0.13.1
rootdir: /private/tmp/nix-build-python3.7-accupy-0.3.2.drv-0/accupy-0.3.2
collected 28 items / 4 deselected / 24 selected

test/test_dot.py .F.........                                             [ 45%]
test/test_sums.py .............                                          [100%]

=================================== FAILURES ===================================
________________________ test_kdot2[1000000000000000.0] ________________________

cond = 1000000000000000.0

    @pytest.mark.parametrize("cond", [1.0, 1.0e15])
    def test_kdot2(cond):
        x, y, ref, _ = accupy.generate_ill_conditioned_dot_product(100, cond)
>       assert abs(accupy.kdot(x, y, K=2) - ref) < 1.0e-15 * abs(ref)
E       AssertionError: assert mpf('0.000000000000000114084570757375210839671410496456786853711854747527032483134235008037649095058441162109375') < (1e-15 * mpf('0.03302009849489094927194480185983845556774002768261744316314525247296751686576499196235090494155883789063'))
E        +  where mpf('0.000000000000000114084570757375210839671410496456786853711854747527032483134235008037649095058441162109375') = abs((0.03302009849489106 - mpf('0.03302009849489094927194480185983845556774002768261744316314525247296751686576499196235090494155883789063')))
E        +    where 0.03302009849489106 = <function kdot at 0x11eede8c0>(array([ 5.10653219e+07, -1.35558140e+04,  1.16214778e+03, -7.64036713e+05,\n       -4.15362782e+02, -5.10518025e+03, -7...8710e+01, -4.28172558e+00,  1.95177498e+00,\n        3.74687011e-01, -2.86712378e+00,  1.99264877e-02,  1.83323108e+00]), array([ 2.28091352e+07, -1.13769239e+04,  9.42159060e+02,  3.52913613e+05,\n        7.06341353e+01, -4.42719693e+03,  2...7027e-01,  2.37823667e+00,  3.99193374e+00,\n       -3.83499060e+00,  1.10856098e+00, -1.75602080e+00,  7.20764693e-01]), K=2)
E        +      where <function kdot at 0x11eede8c0> = accupy.kdot
E        +  and   mpf('0.03302009849489094927194480185983845556774002768261744316314525247296751686576499196235090494155883789063') = abs(mpf('0.03302009849489094927194480185983845556774002768261744316314525247296751686576499196235090494155883789063'))

test/test_dot.py:12: AssertionError
=========================== short test summary info ============================
FAILED test/test_dot.py::test_kdot2[1000000000000000.0] - AssertionError: ass...
================= 1 failed, 23 passed, 4 deselected in 14.64s ==================

A similar build with the same dependencies passes this test on x86_64 linux. I'm currently working on the Nix package for accupy and obviously don't want to ship a broken package. Is this a showstopper?

Relevant dependencies:

/nix/store/9krlzvny65gdc8s7kpb6lkx8cd02c25b-default-builder.sh
/nix/store/x8sbijw81kcc2dcfycy2nm9yxnl5wrma-bash-4.4-p23.drv
/nix/store/87n9yqbs28rr1js0yxfwszdp8ldn0r2j-stdenv-darwin.drv
/nix/store/0nn6gan1407fc1mkh17zcpl1zn6rdxp9-python3-3.7.9.drv
/nix/store/13bk41y1aagcifilqwnr6x6rz0w0b2rf-python-recompile-bytecode-hook.drv
/nix/store/1l0n6pmfgfna1ffgk3g9w0v0shkzzqwv-accupy-0.3.2.tar.gz.drv
/nix/store/6qx7xgsyi25fwyxjcj2dz90k3xxcfs1b-python-imports-check-hook.sh.drv
/nix/store/8bzixyfrrdd3fni0dmh4nw5kknbi5xng-python-remove-bin-bytecode-hook.drv
/nix/store/bk2p76spz2ba25m3iw978m1zpwb6wgmx-hook.drv
/nix/store/g03kwd95zf3wrdp94a1nlvkbav9nlg4j-python-namespaces-hook.sh.drv
/nix/store/k4a197bybzv6jgkk7gjp59fgg64yv1gn-hook.drv
/nix/store/yrkgh4waq5pnv8gpsscfba3r1q0nhl6l-python-remove-tests-dir-hook.drv
/nix/store/bpyv1qgi6k3fjzq7ay4b8n5fp49bfp0h-python3.7-setuptools-47.3.1.drv
/nix/store/32ymxb64r5d31jijzq8h078zgw8pbfag-setuptools-check-hook.drv
/nix/store/4b61dddsswbmzr8xl1qvn4pk2zfmi2g5-python-catch-conflicts-hook.drv
/nix/store/ffx3x2c6chq7164bxsa3p0119mx0mxam-pip-install-hook.drv
/nix/store/mdr3rdy79bj3zvmnrz4k0mm7ipk4skzq-setuptools-setup-hook.drv
/nix/store/ahgmj7j0d95rjbgaq9wrgnv41dxcgb5p-python3.7-numpy-1.19.1.drv
/nix/store/m53dbcxp4mh1jla9q5ymdsvp70zfj8il-eigen-3.3.7.drv
/nix/store/kmx7gnv8nwzpf5wrr730m2jlr4qxvsnj-python3.7-pybind11-2.5.0.drv
/nix/store/4yy64nasiak07w39ghm0scp2x1qyxz7j-python3.7-pyfma-0.1.1.drv
/nix/store/7fpcagcz29gvsi2k4j6dqac6251gxhla-python3.7-matplotlib-3.3.1.drv
/nix/store/jzq99nkn3msjxb30lp1li6dn4wsrahvc-python3.7-mpmath-1.1.0.drv
/nix/store/pcn6qi5z4mhi2sc3mpzk9mlinn927svp-pytest-check-hook.drv
/nix/store/vj9md5s8d8h9w445qgf0syaljkir8rna-python3.7-pipdate-0.5.2.drv

readme images

images/plots in the readme don't seem to be loading

Numpy accuracy

I'm not seeing numpy.sum being that inaccurate, for example:

>>> import numpy as np
>>> import math
>>> np.random.seed(42)
>>> dist = np.hstack([
...     np.random.normal(loc=1, scale=2., size=500_000),
...     np.random.normal(loc=1, scale=.5, size=500_000)
... ])

 # fsum vs. naive
>>> math.fsum(dist) - sum(dist)
5.820766091346741e-10

 # fsum vs. np.sum
>>> math.fsum(dist) - np.sum(dist)
0.0

# Check with 32 bit, fsum vs. naive
>>> math.fsum(dist.astype(np.float32)) - sum(dist.astype(np.float32), np.float32(0))
-20.27040922746528

# Check with 32 bit, fsum vs. np.sum
>>> math.fsum(dist.astype(np.float32)) - np.sum(dist.astype(np.float32))
-0.1454092274652794

# Check with 32 bit fsum vs. 64 bit fsum
>>> math.fsum(dist) - math.fsum(dist.astype(np.float32))
-2.3222528398036957e-06

Failed to compile on Mac OS Mojave 10.14.5 (18F132)

I tried installing accupy and compilation failed:

~/partitioning: python3 -m pip install accupy
Collecting accupy
...
    Complete output from command /usr/local/bin/python3 -u -c "import setuptools, tokenize;__file__='/private/var/folders/_j/_5vgryt961bg8tw_2jvgmc6w0000gp/T/pip-install-t2i_vxvy/pyfma/setup.py';f=getattr(tokenize, 'open', open)(__file__);code=f.read().replace([87/1791]
\n');f.close();exec(compile(code, __file__, 'exec'))" install --record /private/var/folders/_j/_5vgryt961bg8tw_2jvgmc6w0000gp/T/pip-record-38rgrthc/install-record.txt --single-version-externally-managed --compile:
    running install
    running build
    running build_py
    creating build
    creating build/lib.macosx-10.9-x86_64-3.7
    creating build/lib.macosx-10.9-x86_64-3.7/pyfma
    copying pyfma/__init__.py -> build/lib.macosx-10.9-x86_64-3.7/pyfma
    copying pyfma/__about__.py -> build/lib.macosx-10.9-x86_64-3.7/pyfma
    running build_ext
    building '_pyfma' extension
    creating build/temp.macosx-10.9-x86_64-3.7
    creating build/temp.macosx-10.9-x86_64-3.7/src
    gcc -Wno-unused-result -Wsign-compare -Wunreachable-code -fno-common -dynamic -DNDEBUG -g -fwrapv -O3 -Wall -arch x86_64 -g -I/Library/Frameworks/Python.framework/Versions/3.7/include/python3.7m -I/Users/smukherjee/Library/Python/3.7/include/python3.7m -I/Library/Fr
ameworks/Python.framework/Versions/3.7/include/python3.7m -c src/pybind11.cpp -o build/temp.macosx-10.9-x86_64-3.7/src/pybind11.o
    In file included from src/pybind11.cpp:3:
    In file included from /Library/Frameworks/Python.framework/Versions/3.7/include/python3.7m/pybind11/pybind11.h:44:
    In file included from /Library/Frameworks/Python.framework/Versions/3.7/include/python3.7m/pybind11/attr.h:13:
    In file included from /Library/Frameworks/Python.framework/Versions/3.7/include/python3.7m/pybind11/cast.h:13:
    In file included from /Library/Frameworks/Python.framework/Versions/3.7/include/python3.7m/pybind11/pytypes.h:12:
    /Library/Frameworks/Python.framework/Versions/3.7/include/python3.7m/pybind11/detail/common.h:298:17: warning: alias declarations are a C++11 extension [-Wc++11-extensions]
    using ssize_t = Py_ssize_t;
                    ^
    /Library/Frameworks/Python.framework/Versions/3.7/include/python3.7m/pybind11/detail/common.h:299:17: warning: alias declarations are a C++11 extension [-Wc++11-extensions]
    using size_t  = std::size_t;
                    ^
    /Library/Frameworks/Python.framework/Versions/3.7/include/python3.7m/pybind11/detail/common.h:302:6: warning: scoped enumerations are a C++11 extension [-Wc++11-extensions]
    enum class return_value_policy : uint8_t {
         ^
    /Library/Frameworks/Python.framework/Versions/3.7/include/python3.7m/pybind11/detail/common.h:302:12: error: ISO C++ forbids forward references to 'enum' types
    enum class return_value_policy : uint8_t {
               ^
    /Library/Frameworks/Python.framework/Versions/3.7/include/python3.7m/pybind11/detail/common.h:302:32: error: expected unqualified-id
    enum class return_value_policy : uint8_t {
                                   ^
    /Library/Frameworks/Python.framework/Versions/3.7/include/python3.7m/pybind11/detail/common.h:355:15: error: unknown type name 'constexpr'
    inline static constexpr int log2(size_t n, int k = 0) { return (n <= 1) ? k : log2(n >> 1, k + 1); }
                  ^
    /Library/Frameworks/Python.framework/Versions/3.7/include/python3.7m/pybind11/detail/common.h:355:25: error: expected unqualified-id
    inline static constexpr int log2(size_t n, int k = 0) { return (n <= 1) ? k : log2(n >> 1, k + 1); }
                            ^
    /Library/Frameworks/Python.framework/Versions/3.7/include/python3.7m/pybind11/detail/common.h:358:15: error: unknown type name 'constexpr'
    inline static constexpr size_t size_in_ptrs(size_t s) { return 1 + ((s - 1) >> log2(sizeof(void *))); }
                  ^
    /Library/Frameworks/Python.framework/Versions/3.7/include/python3.7m/pybind11/detail/common.h:358:1: warning: inline variables are a C++17 extension [-Wc++17-extensions]
    inline static constexpr size_t size_in_ptrs(size_t s) { return 1 + ((s - 1) >> log2(sizeof(void *))); }
    ^
    /Library/Frameworks/Python.framework/Versions/3.7/include/python3.7m/pybind11/detail/common.h:358:31: error: expected ';' after top level declarator
    inline static constexpr size_t size_in_ptrs(size_t s) { return 1 + ((s - 1) >> log2(sizeof(void *))); }
                                  ^
                                  ;
    /Library/Frameworks/Python.framework/Versions/3.7/include/python3.7m/pybind11/detail/common.h:386:39: error: use of undeclared identifier 'instance_simple_holder_in_ptrs'
            void *simple_value_holder[1 + instance_simple_holder_in_ptrs()];
                                          ^
    /Library/Frameworks/Python.framework/Versions/3.7/include/python3.7m/pybind11/detail/common.h:433:49: error: unknown type name 'type_info'; did you mean 'std::type_info'?
        value_and_holder get_value_and_holder(const type_info *find_type = nullptr, bool throw_if_missing = true);
                                                    ^~~~~~~~~
                                                    std::type_info
    /Library/Developer/CommandLineTools/usr/include/c++/v1/typeinfo:87:29: note: 'std::type_info' declared here                                                                                                                                                      [28/1791]
    class _LIBCPP_EXCEPTION_ABI type_info
                                ^
    In file included from src/pybind11.cpp:3:
    In file included from /Library/Frameworks/Python.framework/Versions/3.7/include/python3.7m/pybind11/pybind11.h:44:
    In file included from /Library/Frameworks/Python.framework/Versions/3.7/include/python3.7m/pybind11/attr.h:13:
    In file included from /Library/Frameworks/Python.framework/Versions/3.7/include/python3.7m/pybind11/cast.h:13:
    In file included from /Library/Frameworks/Python.framework/Versions/3.7/include/python3.7m/pybind11/pytypes.h:12:
    /Library/Frameworks/Python.framework/Versions/3.7/include/python3.7m/pybind11/detail/common.h:436:12: error: unknown type name 'constexpr'
        static constexpr uint8_t status_holder_constructed  = 1;
               ^
    /Library/Frameworks/Python.framework/Versions/3.7/include/python3.7m/pybind11/detail/common.h:436:29: error: expected ';' at end of declaration list
        static constexpr uint8_t status_holder_constructed  = 1;
                                ^
                                ;
    /Library/Frameworks/Python.framework/Versions/3.7/include/python3.7m/pybind11/detail/common.h:437:12: error: unknown type name 'constexpr'
        static constexpr uint8_t status_instance_registered = 2;
               ^
    /Library/Frameworks/Python.framework/Versions/3.7/include/python3.7m/pybind11/detail/common.h:437:29: error: expected ';' at end of declaration list
        static constexpr uint8_t status_instance_registered = 2;
                                ^
                                ;
    /Library/Frameworks/Python.framework/Versions/3.7/include/python3.7m/pybind11/detail/common.h:449:58: warning: alias declarations are a C++11 extension [-Wc++11-extensions]
    template <bool B, typename T = void> using enable_if_t = typename std::enable_if<B, T>::type;
                                                             ^
    /Library/Frameworks/Python.framework/Versions/3.7/include/python3.7m/pybind11/detail/common.h:450:65: warning: alias declarations are a C++11 extension [-Wc++11-extensions]
    template <bool B, typename T, typename F> using conditional_t = typename std::conditional<B, T, F>::type;
                                                                    ^
    /Library/Frameworks/Python.framework/Versions/3.7/include/python3.7m/pybind11/detail/common.h:451:43: warning: alias declarations are a C++11 extension [-Wc++11-extensions]
    template <typename T> using remove_cv_t = typename std::remove_cv<T>::type;
                                              ^
    /Library/Frameworks/Python.framework/Versions/3.7/include/python3.7m/pybind11/detail/common.h:452:50: warning: alias declarations are a C++11 extension [-Wc++11-extensions]
    template <typename T> using remove_reference_t = typename std::remove_reference<T>::type;
                                                     ^
    /Library/Frameworks/Python.framework/Versions/3.7/include/python3.7m/pybind11/detail/common.h:460:10: error: unknown type name 'size_t'
    template<size_t ...> struct index_sequence  { };
             ^
    /Library/Frameworks/Python.framework/Versions/3.7/include/python3.7m/pybind11/detail/common.h:460:17: warning: variadic templates are a C++11 extension [-Wc++11-extensions]
    template<size_t ...> struct index_sequence  { };
                    ^
    /Library/Frameworks/Python.framework/Versions/3.7/include/python3.7m/pybind11/detail/common.h:461:10: error: unknown type name 'size_t'
    template<size_t N, size_t ...S> struct make_index_sequence_impl : make_index_sequence_impl <N - 1, N - 1, S...> { };
             ^
    /Library/Frameworks/Python.framework/Versions/3.7/include/python3.7m/pybind11/detail/common.h:461:20: error: unknown type name 'size_t'
    template<size_t N, size_t ...S> struct make_index_sequence_impl : make_index_sequence_impl <N - 1, N - 1, S...> { };
                       ^
    /Library/Frameworks/Python.framework/Versions/3.7/include/python3.7m/pybind11/detail/common.h:461:27: warning: variadic templates are a C++11 extension [-Wc++11-extensions]
    template<size_t N, size_t ...S> struct make_index_sequence_impl : make_index_sequence_impl <N - 1, N - 1, S...> { };
                              ^
    /Library/Frameworks/Python.framework/Versions/3.7/include/python3.7m/pybind11/detail/common.h:462:10: error: unknown type name 'size_t'
    template<size_t ...S> struct make_index_sequence_impl <0, S...> { typedef index_sequence<S...> type; };
             ^
    /Library/Frameworks/Python.framework/Versions/3.7/include/python3.7m/pybind11/detail/common.h:462:17: warning: variadic templates are a C++11 extension [-Wc++11-extensions]
    template<size_t ...S> struct make_index_sequence_impl <0, S...> { typedef index_sequence<S...> type; };
                    ^
    /Library/Frameworks/Python.framework/Versions/3.7/include/python3.7m/pybind11/detail/common.h:463:10: error: unknown type name 'size_t'
    template<size_t N> using make_index_sequence = typename make_index_sequence_impl<N>::type;
             ^
    /Library/Frameworks/Python.framework/Versions/3.7/include/python3.7m/pybind11/detail/common.h:463:48: warning: alias declarations are a C++11 extension [-Wc++11-extensions]
    template<size_t N> using make_index_sequence = typename make_index_sequence_impl<N>::type;
                                                   ^
    /Library/Frameworks/Python.framework/Versions/3.7/include/python3.7m/pybind11/detail/common.h:467:32: error: expected template parameter
    template <typename ISeq, size_t, bool...> struct select_indices_impl { using type = ISeq; };
                                   ^
    /Library/Frameworks/Python.framework/Versions/3.7/include/python3.7m/pybind11/detail/common.h:467:38: warning: variadic templates are a C++11 extension [-Wc++11-extensions]
    template <typename ISeq, size_t, bool...> struct select_indices_impl { using type = ISeq; };
                                         ^
    /Library/Frameworks/Python.framework/Versions/3.7/include/python3.7m/pybind11/detail/common.h:467:85: warning: alias declarations are a C++11 extension [-Wc++11-extensions]
    template <typename ISeq, size_t, bool...> struct select_indices_impl { using type = ISeq; };
                                                                                        ^
    /Library/Frameworks/Python.framework/Versions/3.7/include/python3.7m/pybind11/detail/common.h:468:11: error: unknown type name 'size_t'
    template <size_t... IPrev, size_t I, bool B, bool... Bs> struct select_indices_impl<index_sequence<IPrev...>, I, B, Bs...>
              ^
    /Library/Frameworks/Python.framework/Versions/3.7/include/python3.7m/pybind11/detail/common.h:468:17: warning: variadic templates are a C++11 extension [-Wc++11-extensions]
    template <size_t... IPrev, size_t I, bool B, bool... Bs> struct select_indices_impl<index_sequence<IPrev...>, I, B, Bs...>
                    ^
    fatal error: too many errors emitted, stopping now [-ferror-limit=]
    15 warnings and 20 errors generated.
    error: command 'gcc' failed with exit status 1

It looks like one of the problems is that my compiler is not expecting C++11. I'm going to try to see if if I can pass a flag to gcc to get it to accept C++ 11.

Recommend Projects

React

A declarative, efficient, and flexible JavaScript library for building user interfaces.
Vue.js

🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
Typescript

TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
TensorFlow

An Open Source Machine Learning Framework for Everyone
Django

The Web framework for perfectionists with deadlines.
Laravel

A PHP framework for web artisans
D3

Bring data to life with SVG, Canvas and HTML. 📊📈🎉

Recommend Topics

javascript

JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
web

Some thing interesting about web. New door for the world.
server

A server is a program made to process requests and deliver data to clients.
Machine learning

Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Visualization

Some thing interesting about visualization, use data art
Game

Some thing interesting about game, make everyone happy.

Recommend Org

Facebook

We are working to build community through open source technology. NB: members must have two-factor auth.
Microsoft

Open source projects and samples from Microsoft.
Google

Google ❤️ Open Source for everyone.
Alibaba

Alibaba Open Source for everyone
D3

Data-Driven Documents codes.
Tencent

China tencent open source team.