[new concept under development] Python's fsum and others

[9il]

See the discussion #2513.

Example 0:

// Precise summation
import std.numeric.summation;
import std.math, std.algorithm, std.range;
auto r1 = iota(  1, 501 ).map!(a => (-1.0)^^a/a);
auto r2 = iota(501, 1001).map!(a => (-1.0)^^a/a);
Summator!double s1 = 0, s2 = 0;
foreach(e; r1) s1.put(e); 
foreach(e; r2) s2 -= e;  //put(-e)
s1 -= s2;
assert(s1.sum == -0.69264743055982025);

Example 1:

import std.algorithm : map, sum;
auto ar = [1, 1e100, 1, -1e100].map!(a => a*10000);
double r = 20000;
assert(r == ar.sum!(Summation.kbn));
assert(r == ar.sum!(Summation.kb2));

Example 2:

// Works like Python's fsum, 
// but  current implementation re-establish special
// value semantics across iterations (i.e. handling -Inf + Inf).

import std.algorithm : map, sum;
import std.range;
import core.stdc.tgmath, std.algorithm, std.range;
auto ar = iota(1000.0).map!(a => 1.7.pow(a+1) - 1.7.pow(a)).chain([-(1.7.pow(1000.0))]);
assert(ar.sum!(Summation.precise)  ==  -1.0);
assert(ar.sum!(real, Summation.precise)  ==  -1.0);

[yebblies]

Does using real for intermediates make it generate incorrect results? Or is it just a performance problem?

[9il]

The results is not exactly incorrect. The value should be rounded to the nearest representable floating-point number using the round-half-to-even rule. So precision loss can happen while truncated result from real to double or float. The main unittest looks like this:

version(X86)
{
    ...
    assert(nextDown(s) <= r && nextUp(s) >= r);
}
else
{
    ...
    assert(s == r);
}

[9il]

Win_32 fails with msg:

Error on line 581: expecting target : dependencies

[yebblies]

Oh ok, makes sense.

Win_32 fails with msg:

You left a merge conflict marker in the makefile.

[9il]

Thanks, fixed.

[9il]

It works!

[kuettler]

Has this been released before? Is an alias needed to maintain compatibility?

[9il]

It will be released in 2.067. An alias is not needed since std.numeric.summation is publicly imported in std.numeric.

[kuettler]

This is amazing and, honestly, I am not qualified to comment on any details. There is a lot in your code to teach me, both about D and numeric. That being said, I am all for including this. However, I'd try to restrict the public import to the numeric package as much as possible, to be able to change implementation details later on.

As far as I can tell, the public functions are sumOfLog2s, geometricMean and of course fsum. Maybe the Summator struct can be private for now? I am guessing here.

On the other hand, if you want to make a large part of the module available for others to use, it could be worthwhile to put this in std.experimental.numeric first.

Again, I really like both the function fsum and its implementation.

[9il]

std.experimental.numeric is a posible solution. I wait for few comments from different reviewers.

The Summator is very useful in realtime apps. I have added following example:

import std.range;
import std.algorithm.mutation : swap;

///Moving mean
class MovingAverage
{
    Summator!double summator;
    double[] circularBuffer;
    size_t frontIndex;

    ///operation without rounding
    void put(double x)
    {
        summator += x;
        swap(circularBuffer[frontIndex++], x);
        summator -= x;
        frontIndex %= circularBuffer.length;
    }

    double avg() @property
    {
        return summator.sum() / circularBuffer.length;
    }

    this(double[] buffer)
    {
        assert(!buffer.empty);
        circularBuffer = buffer;
        summator = 0;
        .put(summator, buffer);
    }
}

/// ma always keeps pricese average of last 1000 elements
auto ma = new MovingAverage(iota(0.0, 1000.0).array);
assert(ma.avg == (1000*999/2) / 1000.0);
/// move by 10 elements
put(ma, iota(1000.0, 1010.0));
assert(ma.avg == (1010*1009/2 - 10*9/2) / 1000.0);

[schuetzm]

Please describe shortly what has to be fixed and why.

[denis-sh]

Nitpick: please, add else \n static assert(0); after static if chains because

• another case could be added in future
• static can easily be missed before if

[dnadlinger]

enum member names are generally lower-cased in Phobos.

[9il]

Fixed

[9il]

Summator can be used to compute precise statistics for data distributed over network. So, a simple representation of internal data are required. Should it be pair

alias Data = Tuple!(ptrdiff_t, "overflowDegree", F, "notFinities", F[], "partialSums");
Data toData(){...}
static typeof(this) createFromData(Data data) {...}

or something else?

[MartinNowak]

What's blocking this?

[schuetzm]

Should this comment be part of the documentation?

[9il]

I think no, it is the same behaviour like for normal ar.map!log2.sum.

[9il]

What's blocking this?

rebased

[9il]

ping @MartinNowak

[MartinNowak]

Can you provide some performance numbers, particularly for the default sum case?
Why do we need so many different algorithms, wouldn't a simple fsum cut it?
Maybe a library would be better suited than adding 1.5KLOC of specialized numeric code to phobos?

[MartinNowak]

You should provide a brief and an extended overview for choosing the most suited algorithms, also explaining their drawbacks. Maybe make a table such as http://dlang.org/phobos/std_algorithm.html.

[9il]

@MartinNowak

1. Yes
2. I will create a table like you suggested. We need Summator for parallel or distributed precise computation. We need precise algorithm because it is really precise. We need Pairswise and Kahan algorithms because they can be used with user defined types like matrixes. Pairswise is faster then Kahan but can be used only with random access ranges. We need KBN because it is faster then Kahan (with fixed fabs for float and double) but cannot be used with user defined types. And KBN2 is just second level of KBN.
3. D definitely needs the simplest way to compute precise sum. So this PR is

• Extension of already included by Andrei Pairwise and Kahan with more precise KBN (and it is faster then Kahan)
• Improved precise Python's fsum implemented as an OutputRange that can be used for precise parallel computations.

The are 1.5KLOC because unittests, special cases for complex numbers, Output Range abstraction, documentation. Andrei is preparing Phobos to work with *.di instead of *.d files. So all unittests and docs would be removed for faster compilation.

[John-Colvin]

Wouldn't it be best if summing could be used like this:

auto walkToLast(R)(R r)
{
    ElementType!R e;
    foreach(el; r) e = el;
    return e;
}
auto mean(R)(R r)
{
    static if (hasLength!R)
    {
        return r.sum / r.length;
    }
    else
    {
        auto lengthAndSum = r.cumulativeSum.enumerate.walkToLast;
        return lengthAndSum[1] / lengthAndSum[0];
    }
}

Background:

• "Reduce" operations are evil, information and structure destroying monsters
• libraries should defer using them until absolutely necessary. Allow them to be deferred to the user of a function/type (e.g. the implementor of mean).
• many/most reductions can, if you break free from traditional functional programming and allow mutable state, be imagined as map operations returning forward (or input if necessary) ranges (e.g. cumulativeSum is a map operation).

Moved example I had to a gist as I don't think it's actually relevant: https://gist.github.com/John-Colvin/8f6b63011e66c38067f5

Here is the important bit

What I would really like to see:

1. Implement all of the summation algorithms (as possible) as output ranges.
2. Add something like this:

/**
 * Takes an input range and an output range with member "state" and 
 * returns a range of the value of state after each element of the input range
 * is put in to the output range. 
 */
struct Cumulative(IR, OR)
if (isOutputRange!OR && hasMember!(OR, "state") && isInputRange!IR)
{
    IR ir;
    OR or;
    private bool empty_ = false;
    auto front() @property
    {
        return or.state;
    }
    bool empty() @property
    {
        return empty_;
    }
    void popFront()
    {
        if(ir.empty) empty_ = true;
        else
        {
            put(or, ir.front);
            ir.popFront();
        }
    }
    static if(isForwardRange!IR && hasMember!(OR, "save"))
        auto save() @property
        {
            auto tmp = this;
            tmp.or = or.save;
            tmp.ir = save;
        }
}
auto cumulative(IR, OR)(IR ir, OR or)
{
    auto tmp = Cumulative!(IR, OR)(ir, or);
    tmp.popFront();
    return tmp;
}

3. Add "save" and "state" to all the Summators, ("state" is just a generic alias for "sum")
4. Add something like this:

template cumulativeSum(F, Summation summation = Summation.appropriate)
{
    F cumulativeSum(Range)(Range r)
    {
        alias Sum = //the relevant output range for the given summation
        return r.cumulative(Sum());
    }
}

[John-Colvin]

No need for sliceable or recursion:

/++
$(LUCKY Pairwise summation) algorithm. Range must be a finite range.
+/
private F sumPairwise(Range, F = Unqual!(ForeachType!Range))(Range data)
if (isInputRange!Range && !isInfinite!Range)
{
    import core.bitop : bsf;
    // Works for r with length < 2^^64, in keeping with the use of size_t
    // elsewhere in std.algorithm and std.range on 64 bit platforms.
    F[64] store = F.max;

    size_t idx = 0;
    foreach (i, el; enumerate(data))
    {
        store[idx] = el;
        foreach (_; 0 .. (i+1).bsf)
        {
            store[idx - 1] += store[idx];
            --idx;
        }
        ++idx;
    }

    F s = store[idx - 1];
    foreach_reverse (i; 0 .. idx - 1)
        s += store[i];

    return s;
}

(https://gist.github.com/John-Colvin/a57d1754f266fba96519)

P.S. there might be opportunities for optimised specialisations here, E.g. reading 32 elements at a time and doing an explicit pairwise sum, or maybe even SIMD.

P.P.S. don't expect the implementation above to give exactly the same result as the recursive implementation, they do slightly different summation trees.

[9il]

Good point, thanks!

[9il]

I will play with @John-Colvin ideas. Lets see how we can improve summation concept.

[John-Colvin]

Note for pairwise sum as an output range

Referring to the eager implementation at https://gist.github.com/John-Colvin/a57d1754f266fba96519

sum/state should return the reverse naïve sum of store[0 .. idx], just like is done in the last stage of the eager implementation.

Note that, because Cumulative doesn't cache a value in popFront, you don't have to pay the cost of this sum for elements that you don't care about the value of. The most important case of this is when calling walkToLast to get the total sum.

[John-Colvin]

Here's a faster pairwise sum: https://gist.github.com/John-Colvin/72a858b1687a2ab49209
It's an order of magnitude better on arrays and simple ranges like a.map!trivialFoo and about 10-20% better on more complicated ones.

Not sure how well it translates to a lazy setting. Having a 16 element cache that you eagerly fill and then exhaust is a possible approach. It's probably not worth doing it though, because in the more complex loops that are typical of lazy, range-based code, the unrolling is likely to be much less valuable.

Nonetheless, it's a nice optimised version for when eager is OK.

[9il]

Could someone please describe this auto-tester error?

: undefined reference to `_D3std7numeric9summation7__arrayZ'
--- errorlevel 1
posix.mak:466: recipe for target 'checkwhitespace' failed

[John-Colvin]

demangled:

checkwhitespace.o(.text.pure nothrow @nogc @safe void std.numeric.summation.Summator!(uint, 0).Summator.put!(uint[]).put(uint[])+0x2f): In function `pure nothrow @nogc @safe void std.numeric.summation.Summator!(uint, 0).Summator.put!(uint[]).put(uint[]):'
: undefined reference to `std.numeric.summation.__array'

[9il]

demangled:

checkwhitespace.o(.text.pure nothrow @nogc @safe void std.numeric.summation.Summator!(uint, 0).Summator.put!(uint[]).put(uint[])+0x2f): In function `pure nothrow @nogc @safe void std.numeric.summation.Summator!(uint, 0).Summator.put!(uint[]).put(uint[]):'
 : undefined reference to `std.numeric.summation.__array'

Thanks @John-Colvin ! @MartinNowak, Does this error message means that checkwhitespace.d was builded with old library *.lib files but with new source code?

[quickfur]

Maybe this is a template emission bug?

[9il]

Maybe this is a template emission bug?

I don't think so.

[9il]

Summation algorithms will be part of std.las.