[python] Numpy: Creating a complex array from 2 real ones?

I want to combine 2 parts of the same array to make a complex array:

Data[:,:,:,0] , Data[:,:,:,1]

These don't work:

x = np.complex(Data[:,:,:,0], Data[:,:,:,1])
x = complex(Data[:,:,:,0], Data[:,:,:,1])

Am I missing something? Does numpy not like performing array functions on complex numbers? Here's the error:

TypeError: only length-1 arrays can be converted to Python scalars

That worked for me:

input:

[complex(a,b) for a,b in zip([1,2,3],[1,2,3])]

output:

[(1+4j), (2+5j), (3+6j)]

If your real and imaginary parts are the slices along the last dimension and your array is contiguous along the last dimension, you can just do

A.view(dtype=np.complex128)

If you are using single precision floats, this would be

A.view(dtype=np.complex64)

Here is a fuller example

import numpy as np
from numpy.random import rand
# Randomly choose real and imaginary parts.
# Treat last axis as the real and imaginary parts.
A = rand(100, 2)
# Cast the array as a complex array
# Note that this will now be a 100x1 array
A_comp = A.view(dtype=np.complex128)
# To get the original array A back from the complex version
A = A.view(dtype=np.float64)

If you want to get rid of the extra dimension that stays around from the casting, you could do something like

A_comp = A.view(dtype=np.complex128)[...,0]

This works because, in memory, a complex number is really just two floating point numbers. The first represents the real part, and the second represents the imaginary part. The view method of the array changes the dtype of the array to reflect that you want to treat two adjacent floating point values as a single complex number and updates the dimension accordingly.

This method does not copy any values in the array or perform any new computations, all it does is create a new array object that views the same block of memory differently. That makes it so that this operation can be performed much faster than anything that involves copying values. It also means that any changes made in the complex-valued array will be reflected in the array with the real and imaginary parts.

It may also be a little trickier to recover the original array if you remove the extra axis that is there immediately after the type cast. Things like A_comp[...,np.newaxis].view(np.float64) do not currently work because, as of this writing, NumPy doesn't detect that the array is still C-contiguous when the new axis is added. See this issue. A_comp.view(np.float64).reshape(A.shape) seems to work in most cases though.

I am python novice so this may not be the most efficient method but, if I understand the intent of the question correctly, steps listed below worked for me.

>>> import numpy as np
>>> Data = np.random.random((100, 100, 1000, 2))
>>> result = np.empty(Data.shape[:-1], dtype=complex)
>>> result.real = Data[...,0]; result.imag = Data[...,1]
>>> print Data[0,0,0,0], Data[0,0,0,1], result[0,0,0]
0.0782889873474 0.156087854837 (0.0782889873474+0.156087854837j)

I use the following method:

import numpy as np

real = np.ones((2, 3))
imag = 2*np.ones((2, 3))

complex = np.vectorize(complex)(real, imag)
# OR
complex = real + 1j*imag

This is what your are looking for:

from numpy import array

a=array([1,2,3])
b=array([4,5,6])

a + 1j*b

->array([ 1.+4.j,  2.+5.j,  3.+6.j])

import numpy as np

n = 51 #number of data points
# Suppose the real and imaginary parts are created independently
real_part = np.random.normal(size=n)
imag_part = np.random.normal(size=n)

# Create a complex array - the imaginary part will be equal to zero
z = np.array(real_part, dtype=complex)
# Now define the imaginary part:
z.imag = imag_part
print(z)

There's of course the rather obvious:

Data[...,0] + 1j * Data[...,1]

If you really want to eke out performance (with big arrays), numexpr can be used, which takes advantage of multiple cores.

Setup:

>>> import numpy as np
>>> Data = np.random.randn(64, 64, 64, 2)
>>> x, y = Data[...,0], Data[...,1]

With numexpr:

>>> import numexpr as ne
>>> %timeit result = ne.evaluate("complex(x, y)")
573 µs ± 21.1 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

Compared to fast numpy method:

>>> %timeit result = np.empty(x.shape, dtype=complex); result.real = x; result.imag = y
1.39 ms ± 5.74 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)