Generic mean function

It seems useful to have a generic function to calculate the mean of all elements in a container, so I wrote one. By default it calculates the arithmetic mean, but should be able to accommodate other means (like geometric) if passed the appropriate functors.

Assumptions:

• value_type can be default-constructed to yield a "zero-value" of its type

I am particularly interested in feedback on:

• The validity of my assumptions
• The mean() interface
• How I validate my iterators
• General code style.

mean.hpp

#ifndef MEAN_HPP
#define MEAN_HPP

#include <algorithm>
#include <functional>
#include <iterator>
#include <numeric>
#include <type_traits>

template <class ForwardIt, class DivOp = std::divides<>, class AddOp = std::plus<>,
class = std::enable_if_t<
std::is_base_of<
std::forward_iterator_tag,
typename std::iterator_traits<ForwardIt>::iterator_category>::value>>
typename std::iterator_traits<ForwardIt>::value_type mean(ForwardIt begin,
ForwardIt end,
DivOp div = DivOp(),
{
using value_type = typename std::iterator_traits<ForwardIt>::value_type;
const auto sum = std::accumulate(begin, end, value_type(), add);
const auto count = std::distance(begin, end);
return div(sum, count);
}

#endif


Below is included some example code to demonstrate the usage of mean(). I'm less interested in feedback on the example code, but it is of course still welcome.

main.cpp

#include "mean.hpp"
#include <iostream>
#include <vector>

struct Point
{
Point() : x(0), y(0) {}
Point(int x, int y) : x(x), y(y) {}

int x, y;
};

Point operator+ (const Point& lhs, const Point& rhs)
{
return Point(lhs.x + rhs.x, lhs.y + rhs.y);
}

int main()
{
// Works on primitives
const std::vector<float> scalars{1, 2, 3, 4, 5, 6, 7, 8, 9};
const auto scalarResult = mean(scalars.begin(), scalars.end());
std::cout << scalarResult << std::endl;

// Works on class types with appropriate operators
const std::vector<Point> points{{1,1}, {2,2}, {3,3}};
const auto pointDivision = [] (const Point& sum, const float count)
{
return Point(sum.x / count, sum.y / count);
};
const auto pointResult = mean(points.begin(), points.end(), pointDivision);
std::cout << pointResult.x << " " << pointResult.y << std::endl;
}


There are a few subtle issues in your code. The naive way of calculating an arithmetic mean by summing everything and dividing by the total number of values can have problems when dealing with a large number of values; overflow becomes very dangerous. Secondly, you're using an enable_if to restrict the iterator category to anything that is not a single-pass iterator. I mention this because both problems can be fixed in the same way: using a better algorithm.

This algorithm comes from Knuth (and originally from Kahan I think). Assuming we add a class ResultType to the list of template arguments, it would look something like:

template <class ForwardIt, class ResultType = double>
ResultType mean(ForwardIt begin, ForwardIt end)
{
ResultType avg{};
int count = 1;
for(auto it = begin; it != end; ++it) {
avg += (*it - avg) / count;
++count;
}
return avg;
}


I've removed the AddOp() and DivOp for simplicity, but putting them back in should be simple enough.

The arithmetic mean is not necessarily a closed operation: the mean of integers is, in general, not an integer. Therefore, you should provide a specialization of the template for integer inputs.