Different shape arrays operations

Question

A bit of background: I want to calculate the array factor of a MxN antenna array, which is given by the following equation:

eq1 Where w_i are the complex weight of the i-th element, (x_i,y_i,z_i) is the position of the i-th element, k is the wave number, theta and phi are the elevation and azimuth respectively, and i ranges from 0 to MxN-1.

In the code I have:

-theta and phi are np.mgrid with shape (200,200) each,

-w_i, and (x,y,z)_i are np.array with shape (NxM,) each

so AF is a np.array with shape (200,200) (sum over i).There is no problem so far, and I can get AF easily doing:

af = zeros([theta.shape[0],phi.shape[0]])
for i in range(self.size[0]*self.size[1]):
            af = af + ( w[i]*e**(-1j*(k * x_pos[i]*sin(theta)*cos(phi) + k * y_pos[i]* sin(theta)*sin(phi)+ k * z_pos[i] * cos(theta))) )

Now, each w_i depends on frequency, so AF too, and now I have w_i with shape (NxM,1000) (I have 1000 samples of each w_i in frequency). I tried to use the above code changing

af = zeros([1000,theta.shape[0],phi.shape[0]])

but I get 'operands could not be broadcast together'. I can solve this by using a for loop through the 1000 values, but it is slow and is a bit ugly. So, what is the correct way to do the summation, or the correct way to properly define w_i and AF ?

Any help would be appreciated. Thanks.

edit The code with the new dimension I'm trying to add is the next:

from numpy import *

class AntennaArray:
    def __init__(self,f,asize=None,tipo=None,dx=None,dy=None):
        self.Lambda = 299792458 / f
        self.k = 2*pi/self.Lambda
        self.size = asize
        self.type = tipo
        self._AF_DATA_SIZE = 200
        self.theta,self.phi =  mgrid[0 : pi : self._AF_DATA_SIZE*1j,0 : 2*pi : self._AF_DATA_SIZE*1j]

        self.element_pos = None
        self.element_amp = None
        self.element_pha = None

        if dx == None:
            self.dx = self.Lambda/2
        else:
            self.dx = dx

        if dy == None:
            self.dy = self.Lambda/2
        else:
            self.dy = dy


        self.generate_array()


    def generate_array(self):


        M = self.size[0]
        N = self.size[1]
        dx = self.dx
        dy = self.dy

        x_pos = arange(0,dx*N,dx)
        y_pos = arange(0,dy*M,dy)
        z_pos = 0

        ele = zeros([N*M,3])

        for i in range(M):
            ele[i*N:(i+1)*N,0] = x_pos[:]


        for i in range(M):
            ele[i*N:(i+1)*N,1] = y_pos[i]


        self.element_pos = ele


        #self.array_factor = self.calculate_array_factor()

    def calculate_array_factor(self):

        theta,phi = self.theta,self.phi

        k = self.k

        x_pos = self.element_pos[:,0]
        y_pos = self.element_pos[:,1]
        z_pos = self.element_pos[:,2]

        w = self.element_amp*exp(1j*self.element_pha)



        if len(self.element_pha.shape) > 1:
            #I have f_size samples of w_i(f)
            f_size = self.element_pha.shape[1]
            af = zeros([f_size,theta.shape[0],phi.shape[0]])
        else:
            #I only have w_i
            af = zeros([theta.shape[0],phi.shape[0]])


        for i in range(self.size[0]*self.size[1]):
        **strong text**#This for loop does the summation over i
            af = af + ( w[i]*e**(-1j*(k * x_pos[i]*sin(theta)*cos(phi) + k * y_pos[i]* sin(theta)*sin(phi)+ k * z_pos[i] * cos(theta))) )
        return af

I tried to test it with the next main

from numpy import *
f_points = 10

M = 2
N = 2

a = AntennaArray(5.8e9,[M,N])

a.element_amp = ones([M*N,f_points])
a.element_pha = zeros([M*N,f_points])

af = a.calculate_array_factor()

But I get

ValueError: 'operands could not be broadcast together with shapes (10,) (200,200) '

Note that if I set

a.element_amp = ones([M*N])
a.element_pha = zeros([M*N])

This works well.

Thanks.

Why the `for` loop to implement the equation? You said you have `numpy` arrays, you don't need the for loop. And on what? Seems you have flattened everything. This is totally not needed. What is `self`? Is that inside a method in a class? If you give us more background on your code, is better. — Valentino, Oct 29 '19 at 23:02
I have updated the post with the complete code. The `for` is for the summation only. The others arrays are only 1D because I am summing over them, and the final shape is the (200,200) of the numpy.mgrid. Thanks. — rdm, Oct 30 '19 at 13:48

score 1 · Accepted Answer · edited Jun 20 '20 at 09:12

I had a look at the code, and I think this for loop:

af = zeros([theta.shape[0],phi.shape[0]])
for i in range(self.size[0]*self.size[1]):
            af = af + ( w[i]*e**(-1j*(k * x_pos[i]*sin(theta)*cos(phi) + k * y_pos[i]* sin(theta)*sin(phi)+ k * z_pos[i] * cos(theta))) )

is wrong in many ways. You are mixing dimensions, you cannot loop that way.
And by the way, to make full use of numpy efficiency, never loop over the arrays. It slows down the execution significantly.

I tried to rework that part.

First, I advice you to not use from numpy import *, it's bad practice (see here). Use import numpy as np. I reintroduced the np abbreviation, so you can understand what comes from numpy.

Frequency independent case

This first snippet assumes that w is a 1D array of length 4: I am neglecting the frequency dependency of w, to show you how you can get what you already obtained without the for loop and using instead the power of numpy.

af_points = w[:,np.newaxis,np.newaxis]*np.e**(-1j*
           (k * x_pos[:,np.newaxis,np.newaxis]*np.sin(theta)*np.cos(phi) + 
            k * y_pos[:,np.newaxis,np.newaxis]*np.sin(theta)*np.sin(phi) + 
            k * z_pos[:,np.newaxis,np.newaxis]*np.cos(theta)
           ))
af = np.sum(af_points, axis=0)

I am using numpy broadcasting to obtain a 3D array named af_points, whose shape is (4, 200, 200). To do it, I use np.newaxis to extend the number of axis of an array in order to use broadcasting correctly. More here on np.newaxis.
So, w[:,np.newaxis,np.newaxis] is an array of shape (4, 1, 1). Similarly for x_pos[:,np.newaxis,np.newaxis], y_pos[:,np.newaxis,np.newaxis] and z_pos[:,np.newaxis,np.newaxis]. Since the angles have shape (200, 200), broadcasting can be done, and af_points has shape (4, 200, 200). Finally the sum is done by np.sum, summing over the first axis to obtain a (200, 200) array.

Frequency dependent case

Now w has shape (4, 10), where 10 are the frequency points. The idea is the same, just consider that the frequency is an additional dimension in your numpy arrays: now af_points will be an array of shape (4, 10, 200, 200) where 10 are the f_points you have defined.

To keep it understandable, I've split the calculation:

#exp_point is only the exponent, frequency independent. Will be a (4, 200, 200) array.
exp_points = np.e**(-1j*
            (k * x_pos[:,np.newaxis,np.newaxis]*np.sin(theta)*np.cos(phi) + 
             k * y_pos[:,np.newaxis,np.newaxis]*np.sin(theta)*np.sin(phi) +
             k * z_pos[:,np.newaxis,np.newaxis]*np.cos(theta)
            ))
af_points = w[:,:,np.newaxis,np.newaxis] * exp_points[:,np.newaxis,:,:]
af = np.sum(af_points, axis=0)

And now af has shape (10, 200, 200).

Now I understand what you mean with flattened. I did not understand the concept of broadcasting until your explanation. Thanks a lot. — rdm, Oct 31 '19 at 00:34

Different shape arrays operations

1 Answers1

Frequency independent case

Frequency dependent case