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