Make `MeanNode` a variety of `ReduceNode`

[arcondello]

We're pretty close to being able to make MeanNode a special case of ReduceNode. You can define a ufunc like

template< DType T>
struct Mean : BinaryFunctionMixin<Mean<T>> {
 /// @copydoc Add::result_type
 using result_type = std::conditional<std::integral<T>, double, T>::type;
 // we use the same convention as NumPy

 /// @copydoc Add::reduction_type
 class reduction_type {
 public:
 reduction_type() = delete;

 reduction_type(T value) noexcept : sum_(value), count_(1) {}

 bool operator==(const result_type& rhs) const {
 return static_cast<result_type>(*this) == rhs;
 }
 bool operator==(const reduction_type& rhs) const {
 return sum_ == rhs.sum_ and count_ == rhs.count_;
 }

 explicit operator result_type() const noexcept { return sum_ / count_; }

 private:
 friend Mean;

 result_type sum_; // could use Kahan summation here if we wanted
 ssize_t count_;
 };

 /// @brief Return the average of `lhs` and `rhs`.
 /// @copydetails Add::operator()
 result_type operator()(const DType auto& lhs, const DType auto& rhs) const noexcept {
 return (lhs + rhs) / 2;
 }
 reduction_type operator()(reduction_type lhs, const DType auto& rhs) const noexcept {
 lhs.sum_ += rhs;
 lhs.count_ += 1;
 return lhs;
 }

 /// @brief Revert an average.
 /// @copydetails Add::inverse()
 static std::optional<result_type> inverse(const DType auto& lhs, const DType auto& rhs) noexcept {
 return 2 * lhs - rhs;
 }
 static std::optional<reduction_type> inverse(reduction_type lhs, const DType auto& rhs) noexcept {
 lhs.sum_ -= rhs;
 lhs.count_ -= 1;
 return lhs;
 }

 static ValuesInfo result_bounds(ValuesInfo lhs, ValuesInfo rhs) {
 return ValuesInfo((lhs.max + rhs.max) / 2, (lhs.min + rhs.min) / 2, false);
 }
 static ValuesInfo result_bounds(ValuesInfo bounds, ssize_t) { return bounds; }
 static ValuesInfo result_bounds(ValuesInfo bounds, limit_type) { return bounds; }

 static constexpr bool associative = false;
 static constexpr bool commutative = true;
 static constexpr bool invertible = true;
};

the hitch is around the definition of associative. It works in the reduction case, but not for the bounds calculation. Needs more thought.

[arcondello]

Incidentally, I think the functional_.hpp ufuncs should use the intervals defined in #550 so rather than having separate result_bounds(...) methods, we can just use the operator()(...) methods and define them for intervals.

[arcondello]

I think the same logic applies to SoftmaxNode no, nevermind