Search code examples
pythonarraysmatlabfor-loopstoring-data

Store data from a triple for loop


I have this code that I wrote in MATLAB to store matrices. I was using cells array, but in Python I don't know how to do this.

Would someone know how can I do this?

The MATLAB code is:

S =[0.5 0.7 0.9 1.1]; % distância entre tx e rx[m]
d = 0.2*ones(1,10); 
h = [ 0 0.1 0.2 0.3 0.4]

n_cam = length(d); % numero de camadas
n_alt = length(h); % numero de alturas
n_S = length(S); % numero de receptores
z = zeros(n_cam,n_alt); % profundidade
Rv_h = zeros(n_S,n_alt);
Rv_v = zeros(n_S,n_alt);
Rv = zeros(n_cam,n_alt);
Rh = zeros(n_cam,n_alt);
S_Rv = cell(1,n_S);
S_Rh = cell(1,n_S);
sigma = 0.3*ones(1,n_cam);
sigmaah = zeros(n_S,n_alt);


for i = 1:n_S
    for  j = 1:n_alt
        for k = 1:n_cam
            z(k,j)= (sum(d(1:k))+h(j))/S(i);
            Rv(k,j) = 1/((4*z(k,j)^2+1)^0.5);
            Rh(k,j) = ((4*z(k,j)^2+1)^0.5)-2*z(k,j);
        end 
        Rv_h(i,j) = 1/((4*(h(j)/S(i))^2+1)^0.5);
        Rh_h(i,j)=((4*(h(j)/S(i))^2+1)^0.5)-2*(h(j)/S(i));
    end
    S_Rv(:,i) = {Rv}; % z para cada camada em cada altura, para cada S
    S_Rh(:,i) = {Rh};
end

for i = 1:n_S
    for  j = 1:n_alt
        Rv = cell2mat(S_Rv(1,i));
        Rh = cell2mat(S_Rh(1,i));
        sigma_ah = sigma(1)*(Rh_h(i,j)-Rh(1,j));
        sigma_av = sigma(1)*(Rv_h(i,j)-Rv(1,j));
        for k = 2:(n_cam-1)
           sigma_ah_ant = sigma_ah;
           sigma_av_ant = sigma_av;
           sigma_ah = sigma_ah_ant + sigma(k)*(Rh(k-1,j)-Rh(k,j));
           sigma_av = sigma_av_ant + sigma(k)*(Rv(k-1,j)-Rv(k,j));
        end
        sigmaah (i,j)  = sigma_ah + sigma(end)*Rh(n_cam-1,j)
        sigmaav (i,j)  = sigma_av + sigma(end)*Rv(n_cam-1,j)
    end
end

I was thinking that in Python I could do something like:

n_S = 4
n_alt = 9
n_cam = 6
Rv =[]
for i in range(1,n_S):
    for j in range(1,n_alt):
        for k in range(1,n_cam):
            z[k][j]= (sum(d[0:k])+h[j])/S[i]
            Rv[i][j][k] = 1/((4*z[k,j]**2+1)**0.5)

But it is not working, the error message I get is

list index out of range.


Solution

  • I would suggest taking use of the broadcasting feature

    The loops can then be replaced to clarify the code :

    from numpy import array, ones
    S =array([0.5, 0.7, 0.9, 1.1])
    d = 0.2*ones((10)); 
    h = array([ 0, 0.1, 0.2, 0.3, 0.4])
    
    z = ((d.cumsum()[:, None] + h).ravel() / S[:, None]).reshape((S.size, d.size, h.size))
    Rv = 1 / (4 * z ** 2 + 1)** .5
    # ...etc
    print(Rv[-1])
    

    Which outputs:

        [[0.93979342 0.87789557 0.80873608 0.73994007 0.67572463]
     [0.80873608 0.73994007 0.67572463 0.61782155 0.56652882]
     [0.67572463 0.61782155 0.56652882 0.52145001 0.48191875]
     [0.56652882 0.52145001 0.48191875 0.4472136  0.41665471]
     [0.48191875 0.4472136  0.41665471 0.38963999 0.36565237]
     [0.41665471 0.38963999 0.36565237 0.34425465 0.32507977]
     [0.36565237 0.34425465 0.32507977 0.30782029 0.29221854]
     [0.32507977 0.30782029 0.29221854 0.27805808 0.26515648]
     [0.29221854 0.27805808 0.26515648 0.25335939 0.24253563]
     [0.26515648 0.25335939 0.24253563 0.23257321 0.22337616]]
    

    Which overlaps with the computation in octave/matlab:

    Rv =
    
       0.93979   0.87790   0.80874   0.73994   0.67572
       0.80874   0.73994   0.67572   0.61782   0.56653
       0.67572   0.61782   0.56653   0.52145   0.48192
       0.56653   0.52145   0.48192   0.44721   0.41665
       0.48192   0.44721   0.41665   0.38964   0.36565
       0.41665   0.38964   0.36565   0.34425   0.32508
       0.36565   0.34425   0.32508   0.30782   0.29222
       0.32508   0.30782   0.29222   0.27806   0.26516
       0.29222   0.27806   0.26516   0.25336   0.24254
       0.26516   0.25336   0.24254   0.23257   0.22338   
    

    This reduces the pyramid of doom and is probably faster than the for-loops due to numpy magic. Edit: formatting Edit2: gave check