rewrite _momentX methods in more functional style

[kalmarek]

depending on the width of vector instruction this is up to 4× faster (N=100), the old implementation matches for k = 6;

@devmotion

[andreasnoack]

I think this is a good change. It requires the minimum Julia version to be 1.6 but I think that is okay.

[devmotion]

I think we should use an improved init value if possible - in many cases the result of sum(...) won't be of the same type as zero(m).

[devmotion]

For instance, we could use

Suggested change

	s = sum(x->abs2(x-m), v, init=zero(m))
	init = (zero(eltype(v)) - zero(m))^2
	s = sum(x->(x-m)^2, v; init=init)

[devmotion]

Similarly, I think we need a different init here. Additionally, @inbounds is unsafe here as i might not be an actual index of wv.

[devmotion]

Suggested change

	s = sum(i -> (@inbounds abs2(v[i] - m) * wv[i]), eachindex(v), init=zero(m))
	init = (zero(eltype(v)) - zero(m))^2 * zero(eltype(wv))
	s = sum(i -> (v[i] - m)^2 * wv[i], eachindex(v, wv); init=init)

[kalmarek]

@devmotion if v and wv do not have the same indices eachindex(v, wv) will throw and therefore it seems that @inbounds here should be safe. How about this:

Suggested change

	s = sum(i -> (@inbounds abs2(v[i] - m) * wv[i]), eachindex(v), init=zero(m))
	s = if isempty(v)
	(zero(eltype(v)) - zero(m))^2*zero(eltype(wv))
	elseif eachindex(v) == eachindex(wv)
	sum(i -> (@inbounds (v[i] - m)^2 * wv[i]), eachindex(v, wv))
	else
	sum(i -> ((v_ - m)^2 * wv_) for (v_, wv_) in zip(v, wv))
	end

[devmotion]

Here as well.

[devmotion]

Suggested change

	s = sum(x->(x-m)^3, v, init=zero(m))
	init = (zero(eltype(v)) - zero(m))^3
	s = sum(x->(x-m)^3, v; init=init)

[devmotion]

Same comments as above:

Suggested change

	s = sum(
	i -> (@inbounds (z = (v[i] - m); z * z * z * wv[i])),
	eachindex(v),
	init=zero(m),
	)
	init = (zero(eltype(v)) - zero(m))^3 * zero(eltype(wv))
	s = sum(
	i -> (v[i] - m)^3 * wv[i],
	eachindex(v, wv);
	init=init,
	)

[devmotion]

Suggested change

	s = sum(x-> (z = x-m; abs2(z*z)), v, init=zero(m))
	init = (zero(eltype(v)) - zero(m))^4
	s = sum(x->(x-m)^4, v; init=init)

[devmotion]

Suggested change

	s = sum(
	i -> (@inbounds (z = (v[i] - m); abs2(z * z) * wv[i])),
	eachindex(v),
	init=zero(m),
	)
	init = (zero(eltype(v)) - zero(m))^4 * zero(eltype(wv))
	s = sum(
	i -> (v[i] - m)^4 * wv[i],
	eachindex(v, wv);
	init=init,
	)

[devmotion]

Suggested change

	s = sum(x -> (x - m)^k, v, init=zero(m))
	init = (zero(eltype(v)) - zero(m))^k
	s = sum(x -> (x - m)^k, v; init=init)

[devmotion]

Suggested change

	s = sum(
	i -> (@inbounds (z = (v[i] - m); z^k * wv[i])),
	eachindex(v),
	init=zero(m),
	)
	init = (zero(eltype(v)) - zero(m))^k * zero(eltype(wv))
	s = sum(
	i -> (v[i] - m)^k * wv[i],
	eachindex(v, wv);
	init=init,
	)

[devmotion]

Since it's unspecified whether init is used at all for non-empty collections, an alternative might be to use the suggested init values only for empty collections and use sum without init keyword argument for non-empty collections (as a side effect this would not even require to bump the Julia version it seems).

[kalmarek]

Since it's unspecified whether init is used at all for non-empty collections, an alternative might be to use the suggested init values only for empty collections and use sum without init keyword argument for non-empty collections (as a side effect this would not even require to bump the Julia version it seems).

I'd proceed with this version then. let me know what do you think about the suggestion I made above