use Base.sum instead of handrolled sum with poorly initialized accumulator

[jrevels]

This results in better type stability in cases where the reduction output is not a Float64. I found this while messing about with the code in this Discourse post and a branch of ForwardDiff which adds stack-allocated, low-dimensional gradients.

Note that I found a bunch of places with the same problem throughout the codebase - all of these spots will likely kill performance with non-Float64 number types like Dual.

Setup:

using BenchmarkTools, ForwardDiff, Distributions

f(x1, x2, D) = loglikelihood(Normal(x1, x2), D)
∇f(x1, x2, D) = ForwardDiff.derivative((y1, y2) -> f(y1, y2, D), (x1, x2))

x1, x2 = rand(2);
D = rand(Normal(5.0, 1.0), 10);

Before this PR:

julia> @benchmark ∇f($x1, $x2, $D)
BenchmarkTools.Trial:
  memory estimate:  1008 bytes
  allocs estimate:  32
  --------------
  minimum time:     704.517 ns (0.00% GC)
  median time:      744.134 ns (0.00% GC)
  mean time:        853.260 ns (12.39% GC)
  maximum time:     20.858 μs (92.76% GC)

After this PR:

julia> @benchmark ∇f($x1, $x2, $D)
BenchmarkTools.Trial:
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     323.714 ns (0.00% GC)
  median time:      323.957 ns (0.00% GC)
  mean time:        325.874 ns (0.00% GC)
  maximum time:     1.221 μs (0.00% GC)

[andreasnoack]

It would be great if you could just remove the _ version. It doesn't really do anything. There is a similar method in multivariates which is also uses a literal for initialization of the accumulator so it would be good to have that one fixed as well.

[ararslan]

The problem with this is that x is a scalar, but _logpdf wants x to be an array. I guess that's why it was using a custom sum here, so it could take slices rather than individual elements. Ref JuliaLang/METADATA.jl#8848 (comment)