On Fri, Sep 4, 2015 at 7:43 AM, Rajaditya Mukherjee < rajaditya.mukherjee@gmail.com> wrote:

​Hi Matt,

I am imagining such a library would be greatly useful for numerical scientists and statisticians for precise computing. Can you give us a toy example of how it works ?

I can't post the code without releasing it but I can briefly describe how it works. The general input-range version works by accumulating a mixed number that represents the mean at any given iteration. The mixed number has a whole number part that is of the value_type of the range, and the fractional part is represented as a numerator that is of the iterator's difference_type and a denominator that is of the iterator's difference type. The initial "value" of this state is {0, 0, 0}, which is indeterminate. At any given iteration, you are guaranteed that the fractional part's denominator is exactly equal to the number of elements that have been iterated over, the fractional part's numerator is in the range [0,denominator), and the whole number is the whole number part of the mean. In other words, the fractional part is never reduced. By having a state represented in this manner you are able to communicate back to the user both the exact average as well as the number of elements in the range (since the number of elements is equal to the denominator's value). This guarantee is also useful as we advance the state. The reason the result is 3 values representing a mixed number and not 2 values simply representing a rational number is because a rational number consisting of 2 values of either the difference type or the value type cannot be guaranteed to be able to represent the exact mean without possibility of overflow. With 3 values of appropriate type, you can guarantee that the mean can be represented since the whole number part will always be bounded by the smallest and largest value in the input range (and remember the whole part's type is the value_type of the range). The components of the fractional part use the difference_type of the range since I always represent the fractional part as N/D where D is the number of elements seen so far and N is in the range [0, D). In this manner, the exact mean is guaranteed to able to be represented at any given iteration, and the accumulation function can advance the state in a manner that will never overflow. You use it like any STL algorithm: ////////// auto result = mean(std::begin(range), std::end(range)); ////////// At the moment I do not have an implementation that works piecemeal as is akin to Boost.Accumulators, however since the implementation of the algorithm really is simply a call to std::accumulate, it should be very simple to expose such a form, and I imagine that it would be especially useful. Some further details on how the state is manipulated with each iteration: At any given point, we have the whole number part, the numerator of the fractional part, and the denominator of the fractional part (remember that the denominator is guaranteed to always be the size of the range up to this point). When advancing to the next state, I do an operation that has the same effect as taking a single, rational version of the mixed number N/D, changing this value to N/(D+1), and then adding in V/(D+1), where V is the next value in the range (remember that here D is the current number of elements that we've iterated over). This isn't directly what I do since I'm dealing with a mixed number instead of a rational number, but it's the same overall effect and is done in a way that cannot overflow. The only thing that is really hairy at this point is that there is mixed-mode arithmetic involved (the value_type and the difference_type need to be operands to arithmetic functions). At the moment I handle the mixed-mode arithmetic by explicitly casting, but ultimately it should probably be a customization point of the algorithm, especially if this is to be useful with user-defined, integer-like value types. -- -Matt Calabrese

On 09/04/2015 08:20 PM, Matt Calabrese wrote:

At any given point, we have the whole number part, the numerator of the fractional part, and the denominator of the fractional part (remember that the denominator is guaranteed to always be the size of the range up to this point). When advancing to the next state, I do an operation that has the same effect as taking a single, rational version of the mixed number N/D, changing this value to N/(D+1), and then adding in V/(D+1), where V is the next value in the range (remember that here D is the current number of elements that we've iterated over). This isn't directly what I do since I'm dealing with a mixed number instead of a rational number, but it's the same overall effect and is done in a way that cannot overflow.

Can you provide an example of how it works? For example, if I have an input array with {10, 20, 30}, how are V, N, and D updated in each iteration?

On 09/05/2015 09:20 PM, Matt Calabrese wrote:

Your sample array isn't too interesting in terms of the state changes, so here's a more complex one for the array {11, 17, 3, 31, 27}

State : 0 0/0 Incremented State: 0 0/1

State : 11 0/1 Incremented State: 5 1/2

State : 14 0/2 Incremented State: 9 1/3

State : 10 1/3 Incremented State: 7 3/4

State : 15 2/4 Incremented State: 12 2/5

State : 17 4/5

I still do not understand how you calculate these state values, but looking at them it becomes clear that the average M[i] at iteration i is given by (assuming W[i] is the first entry in the state tuple at iteration i) M[i] = W[i] + (N[i] / D[i]) I am wondering if you can also calculate the running sum from the state variables by? S[i] = W[i] * D[i] + N[i]

Anyway, it seems like there's interest so I'll try to release it. What exactly is the review process for adding algorithms, or is it still somewhat informal?

If they are to be added to Boost.Algorithm, then you need to get the interest of the maintainer. He will decide what level of review is needed, if any.

On Sat, Sep 5, 2015 at 12:20 PM, Matt Calabrese wrote:

Anyway, it seems like there's interest so I'll try to release it. What exactly is the review process for adding algorithms, or is it still somewhat informal?

The algorithm should be able to be adapted to introduce weights. It would work very similarly, only when advancing the state, it would be based on the weight rather than always effectively using (D+1). Similarly, the denominator would no longer always exactly match the number of elements, but instead it would always match the sum of the weights that were encountered. I haven't implemented this weighted variation, but it shouldn't be too much of a change.

Yes, it's pretty informal. You need to convince the library maintainer (me :-)) that it's interesting and that people would use it. I don't think that's going to be too much of a problem. What I would suggest is that you make a pull request, and then post the URL here, and encourage people to comment/review the code that way. -- Marshall P.S. Sorry for the slow response.

On Thu, Sep 17, 2015 at 9:40 AM, Seth wrote:

On 09/17/2015 04:13 PM, Marshall Clow wrote:

...

On Sat, Sep 5, 2015 at 12:20 PM, Matt Calabrese wrote:

The algorithm should be able to be adapted to introduce weights. It would work very similarly, only when advancing the state, it would be based on the weight rather than always effectively using (D+1). Similarly, the denominator would no longer always exactly match the number of elements, but instead it would always match the sum of the weights that were encountered. I haven't implemented this weighted variation, but it shouldn't be too much of a change. I think Boost Accumulator already does these.

E.g.

http://www.boost.org/doc/libs/1_59_0/doc/html/accumulators/user_s_guide.html...

Yeah, it supports weights, if that's specifically what you're referring to, though I thought that when I skimmed through the mean implementation it was not implemented in a way that avoids I.E. overflow. I could have been mistaken, though. Also, the code just got released, so I'll put it up somewhere this weekend if it's not redundant with what Boost.Accumulators already does somewhere. -Matt Calabrese