Implement accumulate and friends #702

tkf · 2019-12-21T06:57:10Z

This PR implements accumulate, cumsum and cumprod for StaticArray.

It const-props with vectors:

julia> function demov()
           A = SA[1, 2]
           cumsum(A)
       end
demov (generic function with 1 method)

julia> @code_typed optimize=true demov()
CodeInfo(
1 ─     return $(QuoteNode([1, 3]))
) => SArray{Tuple{2},Int64,1,2}

Functions are not completely inlined with matrices. Benchmarks show no allocations:

julia> function demom()
           A = SA[1 2; 3 4]
           cumsum(A, dims=Val(1))
       end
demom (generic function with 1 method)

julia> @code_typed optimize=true demom()
CodeInfo(
1 ─ %1 = Base.add_sum::typeof(Base.add_sum)
│   %2 = invoke Base.var"#kw##accumulate"()($(QuoteNode((dims = Val{1}(),)))::NamedTuple{(:dims,),Tuple{Val{1}}}, StaticArrays.accumulate::typeof(accumulate), %1::typeof(Base.add_sum), $(QuoteNode([1 2; 3 4]))::SArray{Tuple{2,2},Int64,2,4})::SArray{Tuple{2,2},Int64,2,4}
└──      return %2
) => SArray{Tuple{2,2},Int64,2,4}

julia> @btime demom()
  7.390 ns (0 allocations: 0 bytes)
2×2 SArray{Tuple{2,2},Int64,2,4} with indices SOneTo(2)×SOneTo(2):
 1  2
 4  6

(edited by @c42f to add: fixes #730)

coveralls · 2019-12-21T07:44:37Z

Coverage increased (+0.2%) to 81.975% when pulling 8db9ad6 on tkf:accumulate into 26d0e05 on JuliaArrays:master.

c42f · 2020-01-07T08:33:34Z

This looks pretty good to me at a quick read, though I should take a closer look. It's impressive that you can get away with this high level of abstraction.

tkf · 2020-01-07T09:49:08Z

I wrote too many foldls last year 😆

c42f

Sorry I took so long to look at this.

I think it's basically good, though I worry just a bit that it's very different from the rest of the code in this package and is asking rather a lot from the compiler ;-) But hey, we should arguably be moving away from using so many generated functions, so perhaps we should just go with it.

src/mapreduce.jl

c42f · 2020-02-14T00:54:57Z

src/mapreduce.jl

+    data = _map(a, CartesianIndices(a)) do _, CI
+        D = _valof(dims)
+        I = Tuple(CI)
+        J = Base.setindex(I, 1, D)


We already import Base.setindex elsewhere in deque.jl so it should be available. Should we just move that import to the main StaticArrays.jl file and use setindex unqualified here?

c42f · 2020-02-14T00:57:57Z

src/mapreduce.jl

+
+struct _InitialValue end
+
+_maybeval(dims::Integer) = Val(Int(dims))


Maybe it's just me, but I can't help reading this as _maybe_eval :-/ Can we change to _maybe_val?

c42f · 2020-02-14T03:55:54Z

src/mapreduce.jl

+    results = _reduce(
+        a,
+        dims,
+        (init = (similar_type(a, Union{}, Size(0))(), init),),


One thing I noticed here is that for length-0 input, the result has eltype of Union{}. I can see why, but in other cases we do try to to have a "sensible" output eltype for convenience, eg SA{Int}[] .+ SA{Int}[] preserves the Int eltype via some hopefully-not-dubious trickery (#664). Is it possible to do something similar here?

One thing I noticed here is that for length-0 input, the result has eltype of Union{}.

I don't think it's possible to solve this in general as the output value depends on the function op. I think another (somewhat) sensible output is to use the eltype of input array. But this is not always the valid answer; e.g., accumulate(//, [1, 2]) (note that Base's accumulate(//, Int[]) magically returns Rational{Int64}[] as it uses promote_op).

As all static arrays (including MArray) are shape-immutable, I don't think it is so harmful to return Union{} element type.

My preference is:

Union{}

promote_op(op, eltype(input), eltype(input))

eltype(input)

I don't think it's possible to solve this in general as the output value depends on the function op.

This is what #664 hacks around for map by using return_type 😬 (but it's no worse than Base which does the same thing for broadcast and collect on empty containers). I've never been quite sure that's the right choice, but if it can't be made to work or is otherwise terrible I'm ok with Union{} :-) Prior to #664 we used Union{} for map and people mostly seemed ok with it.

Ah, I see. I misread the last half of your first comment.

But I think one difficulty of accumulate compared to map is that output "eltype" is kind of ill-defined (or maybe it's better to say length-dependent). It's conceivable to have "type sequence" like this

op(::A, ::A) :: B op(::B, ::A) :: C op(::C, ::A) :: D op(::D, ::A) :: D # finally "stabilized"

But, even though B, C, D are all different types, promote_type(B, C, D) may return type E that is a concrete type. It is also easy to have "oscillatory" output:

osc_op(::Nothing, ::Any) = missing osc_op(::Missing, ::Any) = nothing

So, maybe we should use something like this?

function f(op, acc, x) T = typeof(acc) while true acc = op(acc, x) T = promote_type(T, typeof(acc)) non_existing_value && return T # `non_existing_value` modeling truncation arbitrary-sized input. end end return_type(f, Tuple{typeof(op), typeof(init), eltype(a)})

Impressively, julia can figure out some non-trivial examples like Base._return_type(f, Tuple{typeof(+), Int, Union{Int,Missing}}) and Base._return_type(f, Tuple{typeof(osc_op), Missing, Int}).

I'm not sure if we need to go this far, though. Base is not doing this for accumulate.

You're right, it's confusing to know what the empty case should do here and it's worse than map. With your fix we have

julia> accumulate((a,b)->(a+b)/2, SA[1,2]) 2-element SArray{Tuple{2},Float64,1,2} with indices SOneTo(2): 1.0 1.5 julia> accumulate((a,b)->(a+b)/2, SA[1]) 1-element SArray{Tuple{1},Int64,1,1} with indices SOneTo(1): 1 julia> accumulate((a,b)->(a+b)/2, SA{Int}[]) 0-element SArray{Tuple{0},Int64,1,0} with indices SOneTo(0)

It seems like we must make a somewhat arbitrary choice for the length-0 case with init. You've got

julia> accumulate((a,b)->(a+b)/2, SA{Int}[1], init=0) 1-element SArray{Tuple{1},Float64,1,1} with indices SOneTo(1): 0.5 julia> accumulate((a,b)->(a+b)/2, SA{Int}[], init=0) 0-element SArray{Tuple{0},Float64,1,0} with indices SOneTo(0)

Overall do you feel like this is an improvement on using Union? I do like that cumsum now preserves numeric eltypes and I feel like that could be helpful for users for the same reason that "fixing" map was.

We do have oddities like

julia> cumsum(Int8[1]) 1-element Array{Int64,1}: 1 julia> cumsum(SA{Int8}[1]) 1-element SArray{Tuple{1},Int8,1,1} with indices SOneTo(1): 1 julia> cumsum(SA{Int8}[1,2]) 2-element SArray{Tuple{2},Int64,1,2} with indices SOneTo(2): 1 3

I think that inconsistency between here and Base might be harmless though, and it's certainly not clear how to "fix" it.

Overall do you feel like this is an improvement on using Union?

That's a tough question... 😅 Personally, I think we should be using Union{} eltype everywhere when it's not possible to figure out eltype without relying on the compiler. However, we don't have a good toolset and coding convention for dealing with Union{} eltype ("mutate-or-widen" style, even at the surface API level).

But I think the consistency of the API within a package and with respect to Base matters. As map is already using return_type, overall, I think using return_type seems to be a decent solution here.

src/mapreduce.jl

Co-Authored-By: Chris Foster <chris42f@gmail.com>

This reverts commit 4ca0144.

c42f · 2020-02-17T01:55:17Z

Looks great. I'd still like to see whether we can get everything inlined for the dims case, but as discussed that can be done separately.

Is this good to go from your end?

tkf · 2020-02-17T02:14:26Z

Yeah, I did everything I wanted to add. By the way, I added the inference-based eltype handling but I forgot to comment that. Did you check it?

Re: dims case, I think it'd be better to optimize things in a different round. I think I'll try simplifying _reduce signature at some point. Though maybe using @generated is the option with least hustle...

c42f

By the way, I added the inference-based eltype handling but I forgot to comment that

Yes I noticed that 👍 I had a play around with it and it seems reasonable enough. I think the zero size choice may be somewhat arbitrary for cases other than init=_InitialValue? I guess there might be some weird sense in taking typeof(init)... considering the left-extended sequence init, init+a[1], init+a[1]+a[2], .... However that seems a bit theoretical rather than practical, given the value of init isn't actually used in the empty case.

I'd be happy to go with what's here.

tkf · 2020-02-17T03:15:11Z

src/mapreduce.jl

+@inline Base.cumsum(a::StaticArray; kw...) =
+    accumulate(Base.add_sum, a; init = Base.reduce_empty(Base.add_sum, eltype(a)), kw...)
+@inline Base.cumprod(a::StaticArray; kw...) =
+    accumulate(Base.mul_prod, a; init = Base.reduce_empty(Base.mul_prod, eltype(a)), kw...)


Re: cumsum(SA{Int8}[1]) #702 (comment)

It's easy to make it Base-compatible. So maybe it's worth doing this?

Oh I guess it's not so bad for + given that already has a specialized add_sum. So sure, if it's easy let's make this compatible.

tkf · 2020-02-17T04:19:17Z

I guess there might be some weird sense in taking typeof(init)...

I considered using typeof(init) but it doesn't work if you want to "adjoin" singleton initial value like we are doing. For example:

struct _CustomInitialValue end

function mapaccumulate(f, op, xs)
    accumulate(xs; init = _CustomInitialValue()) do acc, x
        if acc isa _CustomInitialValue
            f(x)
        else
            op(acc, f(x))
        end
    end
end

c42f · 2020-02-18T05:08:23Z

With the use of reduce_first it looks like you've greatly improved consistency with Base for things like cumsum(SA{Int8}[1]). I think that was the last quibble, so do feel free to merge when you're ready.

tkf · 2020-02-18T05:46:02Z

Yeah, I think it's better this way as we are using return_type anyway. (Due to reduce_first, we can't use eltype anymore even when init is not specified. So, we need to use return_type always.)

I tried to merge this but I don't see the merge button. Also, it looks like I can't push to StaticArrays.jl. Maybe I don't have the right permission?

c42f · 2020-02-18T11:46:49Z

Sorry, just a confusion with organization level permissions. Should be fixed now.

tkf · 2020-02-18T22:05:07Z

Thanks! I just merged it. I appreciate that you reviewed this PR very thoroughly!

c42f · 2020-02-19T01:10:34Z

Very nice — thanks for all the thought you put into addressing my quibbles and for some some nice improvement in functionality!

Implement accumulate and friends

3a26ce9

tkf added 2 commits December 21, 2019 00:05

Run tests for accumulate

801d95c

Skip inference tests in Julia 1.1

4985bcf

c42f mentioned this pull request Feb 14, 2020

non-allocating cumsum for SVector #730

Closed

c42f reviewed Feb 14, 2020

View reviewed changes

tkf and others added 5 commits February 13, 2020 22:28

Update src/mapreduce.jl

4ca0144

Co-Authored-By: Chris Foster <chris42f@gmail.com>

Rename: _maybeval -> _maybe_val

df1caa8

Explain how _map is used from _accumulate

df87a52

Revert: (push(ys, y), y)

b665233

This reverts commit 4ca0144.

Comment on why we use vcat

3db20ef

tkf mentioned this pull request Feb 14, 2020

API request: a function to unwrap Val JuliaLang/julia#34758

Open

Use inference to determine element types

6bf5e05

tkf force-pushed the accumulate branch from 8424200 to 6bf5e05 Compare February 15, 2020 10:06

c42f mentioned this pull request Feb 17, 2020

Consider making push/pushfirst/insert widen the type #732

Open

c42f approved these changes Feb 17, 2020

View reviewed changes

Use reduce_empty in cumsum/cumprod for Array-compatibility

5a7963d

tkf commented Feb 17, 2020

View reviewed changes

Use reduce_first instead of reduce_empty

8db9ad6

tkf merged commit 7a78efd into JuliaArrays:master Feb 18, 2020

This was referenced Feb 27, 2020

Discussion: should we use divide-and-conquer reduce for instruction-level parallelism? #748

Open

Define mapfoldl/foldl for static arrays #750

Merged

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Implement accumulate and friends #702

Implement accumulate and friends #702

tkf commented Dec 21, 2019 •

edited by c42f

Loading

coveralls commented Dec 21, 2019 •

edited

Loading

c42f commented Jan 7, 2020

tkf commented Jan 7, 2020

c42f left a comment

c42f Feb 14, 2020

c42f Feb 14, 2020

c42f Feb 14, 2020

tkf Feb 14, 2020

c42f Feb 14, 2020 •

edited

Loading

tkf Feb 14, 2020 •

edited

Loading

c42f Feb 17, 2020

tkf Feb 17, 2020

c42f commented Feb 17, 2020

tkf commented Feb 17, 2020

c42f left a comment

tkf Feb 17, 2020

c42f Feb 17, 2020

tkf commented Feb 17, 2020

c42f commented Feb 18, 2020

tkf commented Feb 18, 2020

c42f commented Feb 18, 2020

tkf commented Feb 18, 2020

c42f commented Feb 19, 2020


		struct _InitialValue end

		_maybeval(dims::Integer) = Val(Int(dims))

Implement accumulate and friends #702

Implement accumulate and friends #702

Conversation

tkf commented Dec 21, 2019 • edited by c42f Loading

coveralls commented Dec 21, 2019 • edited Loading

c42f commented Jan 7, 2020

tkf commented Jan 7, 2020

c42f left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

c42f Feb 14, 2020 • edited Loading

Choose a reason for hiding this comment

tkf Feb 14, 2020 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

c42f commented Feb 17, 2020

tkf commented Feb 17, 2020

c42f left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

tkf commented Feb 17, 2020

c42f commented Feb 18, 2020

tkf commented Feb 18, 2020

c42f commented Feb 18, 2020

tkf commented Feb 18, 2020

c42f commented Feb 19, 2020

tkf commented Dec 21, 2019 •

edited by c42f

Loading

coveralls commented Dec 21, 2019 •

edited

Loading

c42f Feb 14, 2020 •

edited

Loading

tkf Feb 14, 2020 •

edited

Loading