[Python-ideas] Yet another sum function (fractions.sum)

David Mertz mertz at gnosis.cx
Tue Aug 20 20:14:03 CEST 2013


Oscar's improvement to my binsum() idea, using Stefan's excellent point
that we should just use int.__add__ on the numerators, is quite elegant.
 I.e.:

def binsum(iterable):
    bins = defaultdict(int)
    for num in iterable:
        bins[num.denominator] += num.numerator
    return mergesum([F(n, d) for d, n in bins.items()])

Moreover, looking at Oscar's data, it looks like the improved binsum()
ALWAYS beats rationalsum()[*]

[*] OK, technically there are three cases where this isn't true:
d1000/n=100 and d10/n=1000 where it is about one percent slower (although I
think that is in the noise of timing it); and the "pathological" case of
calculating 'e', where no denominators repeat and where rationalsum() beats
everything else by a large margin at the asymptote (or maybe imergesum()
does, but it is because they behave the same here).

So two thoughts:

(1) Use the much simpler binsum(), which *does* accept an iterator, and in
non-pathological cases will also use moderate memory (i.e. N numbers will
fall into approximately log(N) different denominator classes).

(2) Call this specialized sum something like Fraction.sum().  While the
other functions work on, for example, bigint, the builtin sum() is vastly
better there.  I presume a similar result would be true when concrete
optimizations of "Decimal.sum()" are discussed.

The generic case in Python allows for sum()'ing heterogeneous collections
of numbers.  Putting specialized dunder methods into e.g. Fraction.__sum__
is still going to fail in the general case (i.e. fall back to
__builtins__.sum()).  The special case where you KNOW you are summing a
homogenous collection of a certain style of number is special enough that
spelling it Fraction.sum() or math.fsum() or Decimal.sum() is more explicit
and no harder.

I guess on a small point of symmetry, if the spellings above are chosen,
I'd like it also to be true that:

  math.fsum == float.sum



On Tue, Aug 20, 2013 at 9:39 AM, Oscar Benjamin
<oscar.j.benjamin at gmail.com>wrote:

> On 19 August 2013 10:35, Stefan Behnel <stefan_ml at behnel.de> wrote:
> > David Mertz, 16.08.2013 22:06:
> >> My intuition was that one might do better than mergesum() by binning
> >> numbers.  I.e. this implementation:
> >>
> >>   def binsum(seq):
> >>       bins = dict()
> >>       for f in seq:
> >>           bins[f.denominator] = bins.get(f.denominator, 0) + f
> >>       return mergesum(list(bins.values()))
>
> Good thinking David. Actually for some distributions of inputs this
> massively outperforms the other implementations. One common use case I
> have for Fractions is to check the accuracy of floating point
> computation by repeating the computation with Fractions. In this case
> I would be working with Fractions whose denominators are always powers
> of 2 (and usually only a handful of different values). And as Stefan
> says if we're binning on the denominator then we can make it really
> fast by adding the numerators with int.__add__.
>
> > Simply sorting by denominator gives me a visible advantage over the
> above:
> >
> > def sortsum(seq):
> >     def key(f):
> >         return f.denominator if isinstance(f, F) else 1
> >     seq = sorted(seq, key=key)
> >     if len(seq) < 3:
> >         return sum(seq)
> >     return mergesum(seq)
>
> Interesting.
>
> I've developed an iterator version of mergesum:
>
> def imergesum(iterable):
>     stack = []
>     it = iter(iterable)
>     nextpow2 = 1
>     while True:
>         for n2, val in enumerate(it, 1):
>             if n2 == nextpow2:
>                 stack.append(val)
>                 nextpow2 *= 2
>                 break
>             elif n2 % 2:
>                 stack[-1] += val
>             else:
>                 ishift = -1
>                 while not n2 % 2:
>                     val, stack[ishift] = stack[ishift], val
>                     ishift -= 1
>                     n2 >>= 1
>                 stack[ishift] += val
>         else:
>             return mergesum(stack)  # Could just use sum here
>
> It uses log(N) storage so you'd run out of space on an iterator of say
> 2 ** sys.maxint Fractions but that's unlikely to bother anyone.
>
> I combined the different approaches to make a rationalsum() function
> that is scalable and tries to take advantage of the binsum() and
> sortsum() improvements where possible:
>
> def rationalsum(iterable, tablesize=1000):
>     def binsumgen(iterable):
>         iterator = iter(iterable)
>         bins = defaultdict(int)
>         finished = False
>         denom = itemgetter(0)
>         while not finished:
>             finished = True
>             for n, num in zip(range(tablesize), iterator):
>                 bins[num.denominator] += num.numerator
>                 finished = False
>             if len(bins) >= tablesize or finished:
>                 dn = sorted(bins.items(), key=denom)
>                 yield mergesum([F(n, d) for d, n in dn])
>                 bins.clear()
>     return imergesum(binsumgen(iterable))
>
> I tested this on a suite of different rational sequences (full script
> at the bottom):
> 1) float - Fractions made from floats with power of 2 denominators.
> 2) d1000 - the original numbers with denominators uniform in [1, 1000].
> 3) d10 - like above but [1, 10].
> 4) bigint - big ints rather than Fractions (just for comparison).
> 5) e - the sequence 1 + 1 + 1/2 + 1/3! + 1/4! + ... (this is slow).
>
> In different situations the binsum(), sortedsum() and mergesum()
> functions parform differently (in some cases wildly so). The
> rationalsum() function is never the best but always has the same order
> of magnitude as the best. Apart from the bigint case builtins.sum() is
> always among the worst performing.
>
> Here's the benchmark output (Python 3.3 on Windows XP 32-bit - the
> script takes a while):
>
> Running benchmarks (times in microseconds)...
>
>
> --------------
> float
> ---------------
>      n | sum    | merges | imerge | binsum | sortsu | ration
> ------------------------------------------------------------
>     10 | 282    | 269    | 298    | 238    | 284    | 249
>    100 | 2804   | 2717   | 3045   | 766    | 2689   | 780
>   1000 | 29597  | 27468  | 31275  | 2093   | 27151  | 2224
>  10000 | 2.9e+5 | 2.7e+5 | 3.1e+5 | 12908  | 2.7e+5 | 13122
> 100000 | 3.1e+6 | 2.8e+6 | 3.1e+6 | 1.2e+5 | 2.8e+6 | 1.2e+5
>
>
> --------------
> d1000
> ---------------
>      n | sum    | merges | imerge | binsum | sortsu | ration
> ------------------------------------------------------------
>     10 | 214    | 190    | 203    | 318    | 206    | 343
>    100 | 5879   | 2233   | 2578   | 3178   | 2338   | 3197
>   1000 | 2.6e+5 | 24866  | 29015  | 21929  | 23344  | 21901
>  10000 | 6.1e+6 | 2.7e+5 | 3e+5   | 48111  | 1.8e+5 | 48814
> 100000 | 6.6e+7 | 2.7e+6 | 3e+6   | 1.5e+5 | 1.8e+6 | 3.5e+5
>
>
> --------------
> bigint
> ---------------
>      n | sum    | merges | imerge | binsum | sortsu | ration
> ------------------------------------------------------------
>     10 | 1      | 14     | 19     | 19     | 19     | 36
>    100 | 15     | 54     | 140    | 59     | 72     | 81
>   1000 | 158    | 329    | 1294   | 461    | 499    | 550
>  10000 | 1655   | 2957   | 13102  | 4477   | 4497   | 5242
> 100000 | 16539  | 35929  | 1.3e+5 | 44905  | 55101  | 52732
>
>
> --------------
> e
> ---------------
>      n | sum    | merges | imerge | binsum | sortsu | ration
> ------------------------------------------------------------
>     10 | 155    | 156    | 164    | 262    | 166    | 274
>    100 | 9464   | 2777   | 3190   | 9282   | 2893   | 4079
>   1000 | 1.1e+7 | 8.7e+5 | 9.4e+5 | 1.1e+7 | 8.8e+5 | 9e+5
>
>
> --------------
> d10
> ---------------
>      n | sum    | merges | imerge | binsum | sortsu | ration
> ------------------------------------------------------------
>     10 | 142    | 144    | 156    | 139    | 157    | 153
>    100 | 1456   | 1470   | 1822   | 343    | 1541   | 347
>   1000 | 15080  | 14978  | 16955  | 1278   | 15557  | 1241
>  10000 | 1.5e+5 | 1.5e+5 | 1.6e+5 | 9845   | 1.5e+5 | 10343
> 100000 | 1.5e+6 | 1.5e+6 | 1.7e+6 | 96015  | 1.5e+6 | 99856
>
>
> And here's the script:
>
> # tmpsum.py
>
> from __future__ import print_function
>
> def mergesum(seq):
>     while len(seq) > 1:
>         new = [a + b for a, b in zip(seq[:-1:2], seq[1::2])]
>         if len(seq) % 2:
>             new.append(seq[-1])
>         seq = new
>     return seq[0]
>
> def imergesum(iterable):
>     stack = []
>     it = iter(iterable)
>     nextpow2 = 1
>     while True:
>         for n2, val in enumerate(it, 1):
>             if n2 == nextpow2:
>                 stack.append(val)
>                 nextpow2 *= 2
>                 break
>             elif n2 % 2:
>                 stack[-1] += val
>             else:
>                 ishift = -1
>                 while not n2 % 2:
>                     val, stack[ishift] = stack[ishift], val
>                     ishift -= 1
>                     n2 >>= 1
>                 stack[ishift] += val
>         else:
>             return mergesum(stack)
>
> from collections import defaultdict
>
> def binsum(iterable):
>     bins = defaultdict(int)
>     for num in iterable:
>         bins[num.denominator] += num.numerator
>     return mergesum([F(n, d) for d, n in bins.items()])
>
> from operator import attrgetter
>
> def sortsum(seq):
>     seq = sorted(seq, key=attrgetter('denominator'))
>     if len(seq) < 3:
>         return sum(seq)
>     return mergesum(seq)
>
> from operator import itemgetter
>
> def rationalsum(iterable, tablesize=1000):
>     def binsumgen(iterable):
>         iterator = iter(iterable)
>         bins = defaultdict(int)
>         finished = False
>         denom = itemgetter(0)
>         while not finished:
>             finished = True
>             for n, num in zip(range(tablesize), iterator):
>                 bins[num.denominator] += num.numerator
>                 finished = False
>             if len(bins) >= tablesize or finished:
>                 dn = sorted(bins.items(), key=denom)
>                 yield mergesum([F(n, d) for d, n in dn])
>                 bins.clear()
>     return imergesum(binsumgen(iterable))
>
> sumfuncs = sum, mergesum, imergesum, binsum, sortsum, rationalsum
>
> # Just a quick test
> if True:
>     print('testing', end=' ')
>     from fractions import Fraction as F
>     nums = [F(n, 2*n +1) for n in range(3000)]
>     for n in 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 20, 30, 40, 50, 100, 2005:
>         ns = nums[:n]
>         result = sum(ns)
>         for func in sumfuncs:
>             assert func(ns) == result, func.__name__
>         print('.', end='')
>     print(' passed!')
>
> if True:
>     print('generating data...')
>     from random import randint, gauss, seed
>     from fractions import Fraction as F
>     seed(123456789)
>     nmax = 10 ** 5
>     nums_d1000 = [F(randint(-1000, 1000), randint(1, 1000)) for _ in
> range(nmax)]
>     nums_d10 = [F(randint(-10, 10), randint(1, 10)) for _ in range(nmax)]
>     nums_float = [F(*gauss(0, 1).as_integer_ratio()) for _ in range(nmax)]
>     nums_e = [F(1)]
>     for n in range(1, 1000):
>         nums_e.append(nums_e[-1] / n)
>     nums_bigint = [10**65 + randint(1, 100) for n in range(nmax)]
>     nums_all = {
>         'd1000': nums_d1000,
>         'd10': nums_d10,
>         'float': nums_float,
>         'e': nums_e,
>         'bigint': nums_bigint,
>     }
>
> if True:
>     print('\nRunning benchmarks (times in microseconds)...')
>
>     import timeit
>
>     def mytimeit(stmt, setup):
>         n = 10
>         time = lambda : timeit.timeit(stmt, setup, number=n) / n
>         t = time()
>         if t * n < 1e-1:
>             while t * n < 1e-2:
>                 n *= 10
>             if n == 10:
>                 ts = [t, time(), time()]
>             else:
>                 ts = [time(), time(), time()]
>             t = min(ts)
>         return 1e6 * t
>
>     def fmtnum(n):
>         s = str(int(n))
>         if len(s) > 5:
>             s = '%1.2g' % n
>             if s[-4:-1] == 'e+0':
>                 s = s[:-2] + s[-1]
>         return '%-6s' % s
>
>     gapline = '-' * 60
>     fnames = [f.__name__ for f in sumfuncs]
>     header = '     n | ' + ' | '.join('%-6s' % name[:6] for name in fnames)
>     sum = sum
>     setup = 'from __main__ import ' + ', '.join(['nums'] + fnames)
>     setup += '; nums=nums[:%s]'
>
>     for name, nums in nums_all.items():
>         print('\n\n--------------\n%s\n---------------' % name)
>         print(header)
>         print(gapline)
>         for n in 10, 100, 1000, 10000, 100000:
>             if n > len(nums): break
>             times = []
>             for func in sumfuncs:
>                 stmt = '%s(nums)' % (func.__name__,)
>                 times.append(mytimeit(stmt, setup % n))
>             print(('%6s | ' % n), end='')
>             print(' | '.join(fmtnum(t) for t in times))
>
>
> Oscar
> _______________________________________________
> Python-ideas mailing list
> Python-ideas at python.org
> http://mail.python.org/mailman/listinfo/python-ideas
>



-- 
Keeping medicines from the bloodstreams of the sick; food
from the bellies of the hungry; books from the hands of the
uneducated; technology from the underdeveloped; and putting
advocates of freedom in prisons.  Intellectual property is
to the 21st century what the slave trade was to the 16th.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/python-ideas/attachments/20130820/0f0866fd/attachment-0001.html>


More information about the Python-ideas mailing list