Creating a new array with a different element type?

[ChrisRackauckas]

Usually similar returns the same type, but for SArray it does not. Is there no solution for getting "the same array" with a new element type?

if typeof(A) <: StaticArray
  B = similar_type(A,T)
else
  B = similar(A,T)
end

That's the only solution I could come up with, and it requires separate codepaths every time an SArray is encounters and makes it tough to write generic codes. I am hoping there's a method I missed.

[c42f]

The problem with similar is that it assumes the array is mutable - it gives you something uninitialized and assumes the next step of the user will be to initialize it. For immutable arrays this makes no sense, hence the similar_type API, which lets you first get the type, and then supply the elements to construct that type as a tuple. This isn't really different from calling similar(1:10), which instead of returning a UnitRange returns a Vector.

I don't really understand your example code snippet actually: taking the first branch results in a type, but the second branch returns an instance.

Now, that's not to say we've got all this stuff sorted out - there's certainly several API rough edges where we could do better. Perhaps you could describe your use case in a bit more detail?

[c42f]

By the way, we do have

julia> similar(SVector(1,2,3))
3-element MVector{3,Int64}:
 140462817039632
 140462499519488
               0

Which is still a StaticArray, but which is no longer an immutable one.

[ChrisRackauckas]

Perhaps you could describe your use case in a bit more detail?

I'm not doing anything fancy here. I just want to throw an initialized value into a mutable type to update in the future. For a regular AbstractArray, that's an initialized vector that is mutated. For a number, that's a value that then gets updated each iteration. A StaticArray acts like a number but for some reason you can't just copy with a new type (to get units right) a StaticArray like you can a number or mutable vector.

Yes, that example doesn't seem to work. I instead just have create the vector with changed units when it's a StaticArray by doing the mathematical operation (u/t where u is the SArray of unitful quantities and t is some unitful number) since there seems to be no way other generic way to handle this.

Which is still a StaticArray, but which is no longer an immutable one.

That changes the type.

[ChrisRackauckas]

The code that works for this part is

  if typeof(u) <: AbstractArray && !(typeof(u)<:SArray)
    rate_prototype = similar(u/zero(t))
  else
    rate_prototype = u/zero(t)
  end

StaticArrays seem to require a large amount of special handling in ways that other types don't, so I am wondering if I am missing something.

[c42f]

That changes the type

To the extent that similar should return a mutable array, we have no choice about this. I'm not sure what else we can do if you really want to mutate individual elements.

To get away from mutating individual elements you can write algorithms at a higher level so that you operate on the whole set of vector components at once. In this way you can work with SVector just the way you would with numbers - changing the binding rather than the instance when in a loop. This will generate super fast code for StaticArray, but unfortunately not so much for Base.Array. It will, however, work for both.

As a general comment, we have found that a lot of the implementation of array functionality in Base needed to be reimplemented completely for StaticArrays, because the assumption of mutability was so strongly held.

It would be easier to understand what you're trying to do if your example was a bit more complete (ie, I could actully run it). Here's an example of initializing a vector containing unitful quantities, and stripping them off:

julia> using Unitful

julia> velocity = SVector(10.0u"m/s", 20.0u"m/s")
2-element SVector{2,Quantity{Float64, Dimensions:{𝐋 𝐓^-1}, Units:{m s^-1}}}:
 10.0 m s^-1
 20.0 m s^-1

julia> ustrip.(velocity)
2-element SVector{2,Float64}:
 10.0
 20.0

[ChrisRackauckas]

As a general comment, we have found that a lot of the implementation of array functionality in Base needed to be reimplemented completely for StaticArrays, because the assumption of mutability was so strongly held.

Everything for StaticArrays follows the path for Number and not AbstractArrays. If you treat it like a Number, it "just works" in most codes. If you treat it like an AbstractArray, it doesn't work in most codes. That's one reason why I think that it's just classified in the wrong way, since AbstractArray seems to really mean "use mutability" (whereas all algorithms on Number assume no mutability. There are errors where algorithms assume scalars, but adding appropriate .'s around does not change the compiled code for Number but works with SArrays... so it's clear that's the right path to take with every use of SArray that I've encountered)

It would be easier to understand what you're trying to do if your example was a bit more complete (ie, I could actully run it). Here's an example of initializing a vector containing unitful quantities, and stripping them off:

It's just building an AbstractArray without another element-type. Here's a concrete case where it takes the Number approach:

using Unitful
# Number
A = 10.0u"m/s"
t = 1.0u"s"
B = A./t
# Every AbstractArray Except SArray?
A = [10.0u"m/s", 20.0u"m/s"]
t = 1.0u"s"
similar(A,recursive_eltype(A)/typeof(t)) # Recurse `eltype` for array of arrays, in RecursiveArrayTools
# Non-recursive arrays can just use `eltype`
# SArray
velocity = SVector(10.0u"m/s", 20.0u"m/s")
t = 1.0u"s"
B = A./t

[c42f]

The key piece of information missing from your example is what you're going to do with B::Array after this piece of code. Obviously, since you've called similar and ended up with some uninitialized memory, you plan to fill it later somehow. For immutable containers, this is obviously impossible - you only get one chance to insert the elements, and that's at construction time. To lessen the pain, we've made the high level functions map, broadcast etc work properly without allocations. This is why we were forced to invent similar_type - it gets the type so you can call the right constructor with all the elements already computed. To do this without allocation generally requires some code generation.

Generally we try to make Array and StaticArray subtypes behave the same as multidimensional rectangular containers and linear algebra objects. In this sense, they're nothing like Number other than the fact that they're immutable. This case is no different:

julia> using Unitful

julia> A = [10.0u"m/s", 20.0u"m/s"]
2-element Array{Quantity{Float64, Dimensions:{𝐋 𝐓^-1}, Units:{m s^-1}},1}:
 10.0 m s^-1
 20.0 m s^-1

julia> t = 1.0u"s"
1.0 s

julia> B = A ./ t
2-element Array{Quantity{Float64, Dimensions:{𝐋 𝐓^-2}, Units:{m s^-2}},1}:
 10.0 m s^-2
 20.0 m s^-2

julia> AS = @SVector [10.0u"m/s", 20.0u"m/s"]
2-element SVector{2,Quantity{Float64, Dimensions:{𝐋 𝐓^-1}, Units:{m s^-1}}}:
 10.0 m s^-1
 20.0 m s^-1

julia> BS = AS ./ t
2-element SVector{2,Quantity{Float64, Dimensions:{𝐋 𝐓^-2}, Units:{m s^-2}}}:
 10.0 m s^-2
 20.0 m s^-2

Both the syntax and semantics here are exactly the same.

[andyferris]

I was about to jump in and say exactly what @c42f just wrote.

map, broadcast, reduce, etc are your friends when wanting to deal with arbitrary AbstractArray. Using similar generally requires not only mutability on the output, but frequently also knowledge of what the particular array type's implementation is... e.g. what you write after similar for a SparseMatrix would frequently look quite different to Matrix. However, sum and .* and so-on will work very well for Matrix, SpareMatrix, DMatrix (distributed over many computers), AFMatrix (ArrayFire matrix - on GPU), as well as static arrays like SMatrix and MMatrix.

Is there no solution for getting "the same array" with a new element type?

In general, reinterpret is pretty good for this kind of thing. Unfortunately, it only works for structs when they are stored inside a Base.Array. The good news is the compiler is relatively good at removing useless stack copies for things like B = A ./ 1u"s" (Julia codegen will make the change in element type a no-op, while LLVM often will e.g. avoid dividing by 1 if it is safe). Here's proof, on v0.6:

julia> f(x) = x
f (generic function with 1 method)

julia> g(x) = x ./ 1u"s"
g (generic function with 1 method)

julia> @code_native f(SVector(1.0, 2.0))
	.text
Filename: REPL[4]
	pushq	%rbp
	movq	%rsp, %rbp
Source line: 1
	movups	(%rsi), %xmm0
	movups	%xmm0, (%rdi)
	movq	%rdi, %rax
	popq	%rbp
	retq
	nop

julia> @code_native g(SVector(1.0, 2.0))
	.text
Filename: REPL[5]
	pushq	%rbp
	movq	%rsp, %rbp
Source line: 89
	movups	(%rsi), %xmm0
Source line: 1
	movups	%xmm0, (%rdi)
	movq	%rdi, %rax
	popq	%rbp
	retq
	nop

If this were inside a longer function, it would be a no-op (well, ideally - compilers are strange and hard to predict).

[ChrisRackauckas]

A = rand(10000)
using BenchmarkTools
@benchmark similar($A)
function g(A)
  A ./ 2.0
end
@benchmark g($A)

BenchmarkTools.Trial: 
  memory estimate:  78.20 KiB
  allocs estimate:  2
  --------------
  minimum time:     881.150 ns (0.00% GC)
  median time:      6.130 μs (80.42% GC)
  mean time:        5.639 μs (77.07% GC)
  maximum time:     19.682 μs (85.99% GC)
  --------------
  samples:          4416
  evals/sample:     200

BenchmarkTools.Trial: 
  memory estimate:  78.20 KiB
  allocs estimate:  2
  --------------
  minimum time:     6.513 μs (0.00% GC)
  median time:      8.855 μs (0.00% GC)
  mean time:        15.755 μs (37.04% GC)
  maximum time:     908.523 μs (97.34% GC)
  --------------
  samples:          10000
  evals/sample:     4

I guess if it's just an initialization of arrays phase it's okay to take a performance hit like this, but I am still surprised that there's no easy and generic way to just build an array of similar shape with a different element type and will have to resort to actually doing to numerical calculation.

In general, reinterpret is pretty good for this kind of thing.

You mean, reinterpret(T,copy(A))?

e.g. what you write after similar for a SparseMatrix would frequently look quite different to Matrix.

No. The basic similar in this case works for Array, sparse matrices, GPUArrays, all of the DiffEq modeling arrays (RecursiveArrayTools arrays, MultiScaleArrays, DEDataArrays), MArrays, DistributedArrays, and SharedArrays. It doesn't work for ArrayFire.jl because it's missing a method (if it's useful, I'll do a quick PR to fix. But I am not sure how useful ArrayFire is with GPUArrays being a thing, or if ArrayFire works on v0.6 at all with the broadcast changes and no tests for v0.6). It doesn't work with SArrays by choice because of mutability.

Why does similar have to be about mutability?

Generally we try to make Array and StaticArray subtypes behave the same as multidimensional rectangular containers and linear algebra objects. In this sense, they're nothing like Number other than the fact that they're immutable.

But that's essentially where the code differences come in. AbstractArray paths tend pre-allocate to use mutability:

A .= B.*C

Paths on numbers don't:

A = B.*C # Normally don't `.`, but it's not harmful to do so

(that's a simple example which is fixed by the := suggestion, but generally the choice of f vs f! is not)

In DiffEq, the whole stack "just works" on SArrays if you force it down the path written for Number (more generally, "immutable input types"). Given that DiffEq solvers are solving nonlinear equations for roots and solving linear equations internally, I'm very sure that if there was a dispatch in Optim, IterativeSolvers, and NLSolve which worked directly on Numbers (instead of requiring them to be wrapped in a size-1 array), those same codepaths would work directly on SArray if you put .s in the right spots. That's why I say "they act more like a Number": because code "just works" if you treat them like that (and there were extra .s thrown in there for good measure, which normally compile away). And you agree that code doesn't "just work" if you treat them like an AbstractArray. So I am really surprised you don't agree with me on the category of their actions.

I find it really odd that if you put ApproxFun.jl Fun types, arrays of unitful numbers, bigs, symbolic numbers, etc., arrays on the GPU or distributed across multiple nodes, etc., these will all work in generically written codes, but StaticArrays is the outlier that requires you handle it differently (re-write all initializations as math, i.e. A = B./B, and then directly dispatch on SArray in some form to tell it to use a path for non-immutable types). At least, I don't quite understand how to make generic codes work for them without sacrificing the performance of everything else, while every other choice "just works".

[andyferris]

I feel there is a little bit of confusion going on here - e.g. I don't understand what your example was trying to demonstrate. OTOH, you are quite justified in saying that it's hard to write generic code which works optimally for mutable and immutable objects, because it is!

No. The basic similar in this case works for Array, sparse matrices, GPUArrays, all of the DiffEq modeling arrays (RecursiveArrayTools arrays, MultiScaleArrays, DEDataArrays), MArrays, DistributedArrays, and SharedArrays.

Of course - you are completely correct. However, you would be reducing performance considerably if you just did

function f(a::AbstractArray)
    b = similar(a)
    for i = 1:length(a)
        @inbounds b[i] = 2*a[i]
    end
    return b
end 

compared to map(x->2x, a) for most of those types you listed (yes, I simplified, there is indices and linear vs Cartesian indexing to worry about to make this a bit faster, but I don't think these give you parallelism or preserve sparsity).

Why does similar have to be about mutability?

Well, it's essential - this is really important. See the docstring, for example:

help?> similar
search: similar

  similar(array, [element_type=eltype(array)], [dims=size(array)])

  Create an uninitialized mutable array with the given element type and size, based upon the given source array........

If you return a uninitialized array which is similar to the input, but immutable, then what could you possibly do with it? All you get is uninitialized elements which you can never update. Perhaps it gives you a snapshot of RAM or a bunch of zeros, but in general that array would be completely useless.

One way forward might be to provide some way of initializing the output, so out = similar(func, array, T, newsize) where out[i...] = func(i...). When I tried this in the past, I got a lot of overhead from constructing the closure func, but as the compiler gets better this might turn out to be pretty optimal for certain cases (probably not partial reductions, however).

Everything for StaticArrays follows the path for Number and not AbstractArrays. If you treat it like a Number, it "just works" in most codes. If you treat it like an AbstractArray, it doesn't work in most codes. That's one reason why I think that it's just classified in the wrong way, since AbstractArray seems to really mean "use mutability" (whereas all algorithms on Number assume no mutability. There are errors where algorithms assume scalars, but adding appropriate .'s around does not change the compiled code for Number but works with SArrays... so it's clear that's the right path to take with every use of SArray that I've encountered)

In DiffEq, the whole stack "just works" on SArrays if you force it down the path written for Number (more generally, "immutable input types").

I get what you are saying, but in the AbstractArray interface mutability is most certainly optional. Number has quite a different interface. The difference in programming style has to do with mutability. For Number you are using the "functional programming" approach. As you have noted, you can also follow a functional programming approach with mutable arrays, but you may end up with more allocations.

I find it really odd that if you put ApproxFun.jl Fun types, arrays of unitful numbers, bigs, symbolic numbers, etc., arrays on the GPU or distributed across multiple nodes, etc., these will all work in generically written codes, but StaticArrays is the outlier that requires you handle it differently (re-write all initializations as math, i.e. A = B./B, and then directly dispatch on SArray in some form to tell it to use a path for non-immutable types). At least, I don't quite understand how to make generic codes work for them without sacrificing the performance of everything else, while every other choice "just works".

Yes - I agree that this is a problem. Julia does not currently provide any good way of doing non-allocating generic programming over both mutables and immutables. I've opened a couple of discussions on this before, but not much progress has been made because (I think) it is a hard design problem. My observation is that the core devs were hoping to one day get to (a) elide reallocation of arrays automatically whenever it is safe to do so, so that you can use the same generic (and "functional") code paths for numbers, arrays, etc, and (b) allow an efficient and semantically correct way to "mutate an immutable", e.g. via a Ref which is stack allocated, to allow the mutating codepath for stack-allocated things like SArray. Both require robust escape analysis by the compiler, which will be awesome when it lands, but there's no ETA on these.

Philosophically, one could think of using mutability to avoid extra allocation as just an optimization. In an ideal world, I would suggest writing code in a reasonably functional form, and specialize it for mutable types. Of course, the mutable arrays in Julia are much more numerous and useful, so generally this will be done in reverse. Yes, it is frustrating to add another code path just for immutable arrays, and that Number and StaticArray are different enough that they need seperate code paths, and it's doubly frustrating that Base doesn't even provide a trait for immutability to make this a bit easier! At least you can just specialize on StaticArray and assume they are immutable, for now.

So my advice is if you want to get optimal speed in these packages, it seems to me you have little choice but to "suck it up" and write a third code-path for immutable arrays. I am sorry that it's like this; we are well aware of the difficulty, and a variety of people are looking at how to make this easier in the future.

[ChrisRackauckas]

At least you can just specialize on StaticArray and assume they are immutable, for now.

So my advice is if you want to get optimal speed in these packages, it seems to me you have little choice but to "suck it up" and write a third code-path for immutable arrays. I am sorry that it's like this; we are well aware of the difficulty, and a variety of people are looking at how to make this easier in the future.

Yeah, I guess that is the answer for now: assume StaticArray is the only immutable array type (since there is no trait), specialize on them, and wait for a "true solution" involving some kind of stack-allocated small mutables. I was hoping that wasn't the answer, but it is what it is. Thanks for the discussion: that will definitely give me quite a few tools to make sure this works well.

[andyferris]

Cool - let us know if you run into any more trouble.

[c42f]

By the way, these issues and more are covered back in in #32

I agree it's painful to write efficient generic code with immutable vs mutable arrays right now. Perhaps it would help if we provided a trait from this in StaticArrays, in much the same way that base has separate code paths for IndexLinear and IndexCartesian?

[andyferris]

I was thinking of making a micro-package that defines the trait for all types in Base and which other packages could use/extend.

[ChrisRackauckas]

I was thinking of making a micro-package that defines the trait for all types in Base and which other packages could use/extend.

The only issue with this stuff is that it's only useful if it catches on. If you do this, DiffEq will make use of it, but hopefully packages which make types do. I guess it's useful even if just StaticArrays uses the traits. This is definitely something that would need to be in Base down the line though.

[ChrisRackauckas]

I think similar is just the wrong abstraction for this job. I opened an issue in Base to see if there can be something else:

JuliaLang/julia#22218

FWIW, SharedArrays also have this "problem" with similar returning different types, and likely many others. I for some reason assumed that "same type out" was in the informal definition of similar, but instead it seems built around mutability among other concerns?

[c42f]

Yes, similar doesn't preserve the type, even for types in Base (which is why this whole discussion seemed so confusing at the start). This is what I was referring to when I gave the example of UnitRange.