I’m wondering if it’s possible to assign to the real or imaginary part of a complex array without an explicit loop. For example, I thought this might work, but it doesn’t
julia> z = rand(Complex64, 5);
julia> imag.(z) .= 1:5
5-element Array{Float32,1}:
1.0
2.0
3.0
4.0
5.0
julia> z
5-element Array{Complex{Float32},1}:
0.889114+0.0788732im
0.825836+0.846389im
0.945918+0.32109im
0.635547+0.478854im
0.980431+0.22353im
I guess the assignment is to a new vector, not a view of z.
How about z .+= (1:5)*1im?
No, that adds to the existing value which is not, in general, 0.
Correct. I do think this is a case where you need an explicit loop of some sort. That seems fair since this seems like a little bit of an odd use case. Are you sure there is not some better way of initializing for whatever it is you’re trying to do? I’ve never found myself wanting to do something like this as in my experience the typical pattern is initialize and then do arithmetic just like with real numbers.
You can do:
@. z = complex(real(z), newimag)
where newimag is the new imaginary part(s).
There are lots of software systems that don’t have native support for complex data (notably HDF5), so it’s not uncommon to store the real and imaginary parts as separate bands. Also hardware systems (e.g., radio frequency receivers) might have separate physical channels for the real and imaginary parts, so it does make sense to manipulate them independently in some circumstances.
I didn’t know about the @. macro; I’m sure that’ll come in handy. I also didn’t realize Complex types are immutable, so it’s not possible to assign the real and imaginary parts separately. So this does what I need.
Looking at the output of @code_llvm, I’m still irked by this code doing more than required (load real, store real+imag instead of just store imag). Is there a way get a 2xN Float32 view of a 1xN Complex64 vector? That’s basically what I’d do in C or Numpy.
For that you can use reinterpret, though it will give you a 1x2N Float32 view:
julia> z = rand(Complex64, 2);
julia> re_im = reinterpret(Float32, z)
4-element Array{Float32,1}:
0.187461
0.717338
0.352635
0.539329
You can certainly reinterpret an array of Complex as an array of floats:
julia> z = rand(Complex64, 5);
julia> reinterpret(Float32, z, (2, 5))
2×5 Array{Float32,2}:
0.890572 0.2748 0.557549 0.278212 0.133641
0.895787 0.819513 0.132962 0.321889 0.638579
You can even assign to the reinterpreted vector to mutate the underlying data, but I’m not entirely sure it’s a good idea:
julia> reinterpret(Float32, z, (2, 5))[2, :] = 0.5
0.5
julia> z
5-element Array{Complex{Float32},1}:
0.890572+0.5im
0.2748+0.5im
0.557549+0.5im
0.278212+0.5im
0.133641+0.5im
This should be fixed eventually: WIP: Make mutating immutables easier by Keno · Pull Request #21912 · JuliaLang/julia · GitHub
Thanks everyone for the helpful tips. On a related note, I tried comparing the explicit loop to the vector syntax, e.g.,
T = Complex64
m, n = (5000, 5000)
x = rand(T, (m,n))
y = zeros(T, m)
@time for j=1:n
for i=1:m
y[i] = x[i,j]
end
end
@time for j=1:n
y .= x[:,j]
end
and I was surprised the explicit loop is vastly slower on my machine:
$ julia test.jl
3.635253 seconds (94.91 M allocations: 1.787 GiB, 5.77% gc time)
0.149731 seconds (73.28 k allocations: 193.852 MiB, 10.04% gc time)
What’s going on? Am I allocating a new object every time I index a matrix?
You’re timing in global scope. You should wrap your code in a function or, even better, just use GitHub - JuliaCI/BenchmarkTools.jl: A benchmarking framework for the Julia language
Benchmarking with BenchmarkTools shows that your explicit loop is actually faster:
julia> function f1(x, y, n, m)
for j=1:n
for i=1:m
y[i] = x[i,j]
end
end
end
f1 (generic function with 1 method)
julia> function f2(x, y, n)
for j=1:n
y .= x[:,j]
end
end
f2 (generic function with 1 method)
julia> using BenchmarkTools
julia> @benchmark f1($x, $y, $n, $m)
BenchmarkTools.Trial:
memory estimate: 0 bytes
allocs estimate: 0
--------------
minimum time: 20.441 ms (0.00% GC)
median time: 20.733 ms (0.00% GC)
mean time: 20.810 ms (0.00% GC)
maximum time: 22.094 ms (0.00% GC)
--------------
samples: 241
evals/sample: 1
julia> @benchmark f2($x, $y, $n)
BenchmarkTools.Trial:
memory estimate: 191.26 MiB
allocs estimate: 19489
--------------
minimum time: 42.840 ms (7.15% GC)
median time: 43.673 ms (8.05% GC)
mean time: 43.704 ms (7.84% GC)
maximum time: 47.225 ms (6.91% GC)
--------------
samples: 115
evals/sample: 1
That’s because x[:, j] makes a copy of that slice of the array. You can make a view instead, which doesn’t copy the data:
julia> function f3(x, y, n)
for j=1:n
y .= @view(x[:, j])
end
end
f3 (generic function with 1 method)
julia> @benchmark f3($x, $y, $n)
BenchmarkTools.Trial:
memory estimate: 0 bytes
allocs estimate: 0
--------------
minimum time: 13.306 ms (0.00% GC)
median time: 13.952 ms (0.00% GC)
mean time: 14.013 ms (0.00% GC)
maximum time: 15.861 ms (0.00% GC)
--------------
samples: 357
evals/sample: 1
That’s even better than your explicit loop, but you can make your loop just as fast by turning off bounds checks (of course, this is only OK if you actually check the sizes of your arrays to ensure you won’t access out of bounds):
julia> function f4(x, y, n, m)
for j=1:n
for i=1:m
@inbounds y[i] = x[i,j]
end
end
end
f4 (generic function with 1 method)
julia> @benchmark f4($x, $y, $n, $m)
BenchmarkTools.Trial:
memory estimate: 0 bytes
allocs estimate: 0
--------------
minimum time: 13.205 ms (0.00% GC)
median time: 13.421 ms (0.00% GC)
mean time: 13.448 ms (0.00% GC)
maximum time: 14.155 ms (0.00% GC)
--------------
samples: 372
evals/sample: 1
I’m learning a lot here, thanks! I knew about @view but it made no difference in my initial test. If I wrap the test in a function like
function main()
...
end
main()
then I get similar results to @rdeits. What’s special about global scope with respect to speed/memory? For simple programs I wouldn’t have thought it necessary to ensure all number crunching occurs in a function. Or is this an issue with @time specifically?
edit: Oh it turns out I should rtfm. Thanks!
A global variable might have its value, and therefore its type, change at any point. This makes it difficult for the compiler to optimize code using global variables. Variables should be local, or passed as arguments to functions, whenever possible. Any code that is performance critical or being benchmarked should be inside a function.