unrecognized function or variable in parfor

Hello everyone!

I am running the following code and getting this error:

Error using Coverage_chetlur>(parfor supply)

Unrecognized function or variable 'ZZ'.

Error in Coverage_chetlur (line 173)

parfor B = 1:2

Could you please help me fixing it? I am not sure why it appears since I am defining ZZ after having initialized the parfor loop.

PS: Please note that the code works perfectly if I use the for loop instead of the parfor one.

tic
Pi = 3.14159 ;
nIntervals = 3   ;
nPoints = 3 ;

%% PARAMETERS
lambdaT = 0.5e-5  ; 
alpha = [3.5 , 3.5] ; 
m = [2 , 1] ; 
    iterations = 1e0   ;  
    p= [1.5, 10] ;  
    NA = 2   % Better not to exceed NA=20
    Rd  = 5e2 ;
    R_T = Rd + 19e3 ;  
    infty = 1e9 ;   
    zero = 1e-6  ;     
tau = 0.31623 ; 
sigma_n2 = 1e-12 ;    
R_U = linspace(0,1, 10) * Rd ;          R_U = 0.7 * Rd ; 
R_Plus = Rd+R_U-zero ;
R_Minus = Rd-R_U+zero ;
H = 1e2 ;
urbL = 3 ;
eta = [0.9772  0.007943  1 ;  0.7943  0.01  1 ;  0.6918  0.005  1 ;  0.5888  0.0004  1]  ;   
eta = [0.9772    1 ;  0.7943    1 ;  0.6918    1 ;  0.5888    1]  ;   eta(:,2) = eta(:,1) ;  %%% NOT CONSIDERING NLoS
Sa = [4.88 9.6117 12.0810 27.304]' ;   Sb = [0.429 0.1581 0.1139 0.0797]' ;    Sa = Sa(urbL) ;   Sb = Sb(urbL) ;
xi = eta(urbL,:) .* p ;   
D{1,1} = @(z) sqrt(z.^2+H^2) ;
D{1,2} = @(z) (xi(2)/xi(1)).^(1/alpha(2)) .* ( sqrt(z.^2+H^2).^(alpha(1)/alpha(2)) ) ;
D{2,1} = @(z) (xi(1)/xi(2)).^(1/alpha(1)) .* z.^(alpha(2)/alpha(1)) .* (z>D{1,2}(H)) ...
                    + H .* (z<= D{1,2}(H)) ;
D{2,2} = @(z) z ;
    for B = [1,2]
    Z{B,1} = @(z) sqrt(D{B,1}(z).^2-H^2) ;
    Z{B,2} = @(z) D{B,2}(z) ;
    end
 
for u = 1:length(R_U)
Ru = R_U(u) ;   Rplus = Rd+Ru-zero ;   Rminus = Rd-Ru+zero ;
    for it = 1 : iterations                                 
    thetaA = 2*pi*rand(NA,1) ;
    rA  = sqrt( rand(NA,1)*(Rd^2 - 0^2) + 0^2 ) ;   
    coordA = rA .* [cos(thetaA),  sin(thetaA)]   ; 
    Zs{1} = sqrt( (Ru-coordA(:, 1)).^2 + coordA(:, 2).^2 )  ;
    ED{1,it} = sqrt(Zs{1}.^2 + H^2) ;
    EDmin{1}(it) = min(ED{1,it})  ; 
    coord{1} = coordA ;  
    
    nT = poissrnd (lambdaT * pi * (R_T^2-Rd^2)) ;
    thetaT = rand(nT,1) * 2*pi ;
    rT  = sqrt(rand(nT,1)*(R_T^2-Rd^2) + Rd^2) ;  
    coord{2} = [ rT.*cos(thetaT) ,  rT.*sin(thetaT) ] ;
    zT = sqrt( (Ru-coord{2}(:, 1)).^2 + coord{2}(:, 2).^2 )  ;
    ED{2,it} = sqrt( (Ru-coord{2}(:, 1)).^2 + coord{2}(:, 2).^2 ) ;      Zs{2} = ED{2,it} ;
    EDmin{2}(it) = min(ED{2,it})  ;
    end
    

        for B = [1,2]
        mu{B} = @(z) m(B).* tau./xi(B).* D{B,B}(z).^alpha(B) ;
        end
    
%% COVERAGE SIMULATION      
  for it = 1:iterations
       for B = [1,2]
       pTAG(B) = xi(B) * EDmin{B}(it)^-alpha(B) ;  % Conditional received power 
       Hc{B} = gamrnd( m(B),1/m(B),  1,length(ED{B,it}) )   ;
       end
    tag(it) = 1*(pTAG(1)>pTAG(2)) + 2*(pTAG(2)>pTAG(1)) ;   
    receivedPower = Hc{tag(it)}( find(ED{tag(it), it}==EDmin{tag(it)}(it)) ) * xi(tag(it))*EDmin{tag(it)}(it)^-alpha(tag(it)) ;
    interference(it) = sum(Hc{1}'.*xi(1).*ED{1,it}.^-alpha(1)) + sum(Hc{2}'.*xi(2).*ED{2,it}.^-alpha(2)) -  receivedPower ;
    SINR(it) = receivedPower/(interference(it)+sigma_n2) ;
  end
Pc_s(u) =  sum(SINR >= tau) / iterations ;
              
%% ANALYSIS
    Betai = @(z) real( (z>0).*acos( (z.^2 + Ru.^2 - Rd.^2) ./ ( (z>0).*2.*z.*Ru + (z==0)) ) )  ;
    dBetai = @(z) (z>0 ).*(Ru.^2 - z.^2 - Rd.^2) ./ ( (z>0) .* (z .* sqrt(4*z.^2 .* Ru.^2 - (z.^2 + Ru.^2 - Rd.^2).^2)) + (z==0) )  ;
    zX = @(beta) Ru.*cos(beta) + sqrt(Rd.^2 - Ru.^2.*sin(beta).^2) ; 
    rW = @(w,beta) sqrt(Ru^2 + w.^2 - 2*Ru.*w.*cos(beta)) ;
    A0 = Pi* Rd^2 ;
    IND{1} = @(z) z <= Rminus ;
    IND{2} = @(z) (Rminus < z) & (z < Rplus) ;
    IND{3} = @(z) z >= Rplus ;
Sigma_i{1} = @(z) pi*z.^2 ;
Sigma_i{2} = @(z) Betai(z).*z.^2 + integral(@(beta) zX(beta).^2, Betai(z),pi,'ArrayValued',1, 'RelTol',1e-4, 'AbsTol',1e-6) ;
dSigma_i{1} = @(z) 2*pi*z ; 
dSigma_i{2} = @(z) 2*z.*Betai(z) + z.^2 .* dBetai(z) - zX(Betai(z)).^2 .* dBetai(z) ; 
Sigma = @(z) Sigma_i{1}(z).*IND{1}(z) + Sigma_i{2}(z).*IND{2}(z) + (A0+zero) .* IND{3}(z) ;
dSigma = @(z) dSigma_i{1}(z).*IND{1}(z)+dSigma_i{2}(z).*IND{2}(z) ;
F_OmegaA= @(w) Sigma(w) / (pi*Rd^2) ;                             Fbar_OmegaA = @(w) 1-F_OmegaA(w) ;                                f_OmegaA = @(w) dSigma(w) / A0 ;
Fbar_Z{1} = @(z_) arrayfun(@(z) Fbar_OmegaA(z).^NA, z_) ;                                    
Fbar_Z{2} = @(z_) arrayfun(@(z) exp( - lambdaT .* integral(@(beta) integral(@(w) (rW(w,beta)>Rd) .* w,  0, z, 'arrayvalued',1,'RelTol',1e-5, 'AbsTol',1e-8), 0,2*pi, 'arrayvalued',1,'RelTol',1e-5, 'AbsTol',1e-8) ) , z_) ;
F_Z{1} = @(z_) arrayfun(@(z) 1-Fbar_Z{1}(z), z_) ;      
f_Z{1} = @(z_) arrayfun(@(z) NA * Fbar_OmegaA(z).^(NA-1) .* f_OmegaA(z) , z_) ; 
F_Z{2} =@(z_) arrayfun(@(z) 1-Fbar_Z{2}(z), z_) ;
f_Z{2} = @(z_) arrayfun(@(z) lambdaT .* Fbar_Z{2}(z) .* z .* integral(@(beta) (rW(z,beta)>Rd), 0,2*pi,'arrayvalued',1,'RelTol',1e-5, 'AbsTol',1e-8), z_) ;

%% LAPLACE TRANSFORMS
f_hatOmegaA = @(w,z) dSigma(w) ./ (A0-Sigma(z)) ; 
I_{1} = @(s,z) ( m(1) ./ (m(1) + xi(1) * s .* D{1,1}(z).^-alpha(1)) ).^m(1) ;  
I_{2} = @(s,w) 1 - ( m(2) / (m(2)+xi(2) * s .* D{2,2}(w).^-alpha(2)) ).^m(2) ; 
    for B = [1,2]
    hatUps{B,1} = @(s_,z_) arrayfun(@(s,z) integral(@(w) I_{1}(s,w) .* f_hatOmegaA(w,Z{B,1}(z)), Z{B,1}(z), Rplus,'arrayvalued',1, 'RelTol',1e-4, 'AbsTol',1e-8), s_,z_) ;      
    LapBC{B,1} = @(s_,z_) arrayfun(@(s,z) hatUps{B,1}(s,z).^(NA-(B~=2)), s_,z_) ;
    LapBC{B,2} = @(s_,z_) arrayfun(@(s,z) exp( -lambdaT * integral(@(beta) integral(@(w) (rW(w,beta) > Rd) .* I_{2}(s,w) .* w, Z{B,2}(z),infty, ...
                                                     'arrayvalued',1, 'RelTol',1e-4, 'AbsTol',1e-8), 0,2*pi,'arrayvalued',1, 'RelTol',1e-4, 'AbsTol',1e-8) ) , s_,z_)  ; 
    
    LapIII{B,2} = @(s,z) exp(-2*pi*lambdaT * integral(@(x) I_{2}(s,x) .* x,  Z{B,2}(z),infty,'arrayvalued',1) )   ; 
    LapII{B,2} = @(s,z) LapIII{B,2}(s,z) .* exp( lambdaT * integral(@(beta) integral(@(x) I_{2}(s,x) .* x,  Z{B,2}(z),zX(beta),'ArrayValued',1),  -Betai(Z{B,2}(z)),Betai(Z{B,2}(z)),'ArrayValued',1) )  ;  
%     LapI{B,2} = @(s,z) LapIII{B,2}(s,z) .* exp( lambdaT * integral(@(beta) integral(@(x) I_{2}(s,x) .* x,  Z{B,2}(z),zX(beta),'ArrayValued',1),  -pi,pi,'ArrayValued',1) )  ; 
    LapBC{B,2} = @(s,z) LapIII{B,2}(s,z).*(Z{B,2}(z)>=Rplus) + LapII{B,2}(s,z).*( (Rminus<=Z{B,2}(z)) & (Z{B,2}(z)<Rplus) ) ; %+ LapI{B,2}(s,z).*(Z{B,2}(z)<=Rd-Ru(u))  ;
    
    LapB{B} = @(s_,z_) arrayfun(@(s,z) LapBC{B,1}(s,z) .* LapBC{B,2}(s,z) , s_,z_)  ;             
    end

%% ASSOCIATION PROBABILITIES
R_B = [zero,Rplus ; zero,Rplus; Rminus,infty] ;
a{1} = @(z) Fbar_Z{2}(Z{1,2}(z)) ; 
a{2} = @(z) Fbar_Z{1}(Z{2,1}(z)) ;
    for B = [1,2]
    A(B) = integral(@(z) f_Z{B}(z) .* a{B}(z), R_B(B,1), R_B(B,2), 'arrayvalued',1,'RelTol',1e-4, 'AbsTol',1e-8) ;
%% COVERAGE PROBABILITY
    Pcond{B} = @(z_) arrayfun(@(z) 0,  z_) ;
    epsilon(B) = (factorial(m(B)))^(-1/m(B)) ;
  for k = 1:m(B)
        Pcond{B} = @(z_) arrayfun(@(z) Pcond{B}(z) + nchoosek(m(B),k) * (-1)^(k+1) * exp(-k * epsilon(B) * mu{B}(z) * sigma_n2) .* LapB{B}(k * epsilon(B) *mu{B}(z),z),  z_) ; 
  end
    end
%% ANALYTICAL INTEGRATION
Pc_a(u) = 0 ;

 %% MONTECARLO INTEGRATION
 for B = 1:2
    zAxis{B} = logspace( log10(R_B(B,1)), log10(R_B(B,2)), nIntervals) ; 
    prod{B} = zeros(1, nPoints) ; 
    F_zAxis{B} = F_Z{B}(zAxis{B}) ; 
        parfor nP = 1 : nPoints
               nP
        U = rand ;
        diff = abs(U - F_zAxis{B}) ;
        mindiff = min(diff) ;
        index = find(diff==min(diff)) ;
        ZZ{B}(nP) = zAxis{B}(index(1)) ; % We use index(1) because there might be many minimum values        
        prod{B}(nP) = a{B}( ZZ{B}(nP) ) .* Pcond{B}( ZZ{B}(nP) ) ; 
        end
    pc_a{B} = sum(prod{B}) / nPoints ;   
 end
 Pc_a(u) = pc_a{1} + pc_a{2} 
end
toc

[Geoff Hayes]

Maurilio - where is ZZ defined?

[Maurilio Matracia]

Inside Montecarlo integration section, almost at the end

Answer 1

[Raymond Norris]

I don't have much time to look at this closer, so I'll give a couple of quick thoughts

1. For starters, I would suggest going through and addressing the Code Analyzer prompts (the orange markings). There are quite a few cases of preallocation that ought to be done.
2. Without running your code, I can't tell for sure, but I don't think you need to use cell arrays as much as you are.
3. To address the parfor, I would suggest refactoring your code as such

    %% MONTECARLO INTEGRATION
    prod = zeros(2, nPoints);
    for B = 1:2
        zAxis{B} = logspace( log10(R_B(B,1)), log10(R_B(B,2)), nIntervals) ;
        F_zAxis{B} = F_Z{B}(zAxis{B}) ;
        parfor nP = 1 : nPoints
            prod(B,nP) = unit_of_work(F_zAxis,zAxis,a,Pcond,B,nP) ;
        end
        pc_a{B} = sum(prod{B}) / nPoints ;
    end
    Pc_a(u) = pc_a{1} + pc_a{2}

Notice a couple of things

• I'm switching prod from a cell array to a numeric array
• I'm preallocating prod
• I'm putting everything in the parfor into a subfunction (see below)
• I didn't change how B is indexed (e.g. prod{B}), you'll need to do that if you stick with numeric arrays
• I would suggest a different variable name then prod, as this referrs to the built-in function for multiplication of elements of a variable

unit_of_work is defined as follows

function result = unit_of_work(F_zAxis,zAxis,a,Pcond,B,nP)

nP
U = rand ;
diff = abs(U - F_zAxis{B}) ;
index = find(diff==min(diff)) ;
ZZ{B}(nP) = zAxis{B}(index(1)) ; % We use index(1) because there might be many minimum values
result = a{B}( ZZ{B}(nP) ) .* Pcond{B}( ZZ{B}(nP) ) ;

end

You might have to clean something up a bit (i.e. cell array to numeric), but this ought to get you started.

[Maurilio Matracia]

Thank you very much. Actually I just had to initialize ZZ and zAxis to make it work, do you know why?

[Image Analyst]

How would it know what ZZ is otherwise? How would it suppose to know that ZZ is a cell array with at least B cells? And if you're going to assign the nP index to the contents of the B'th cell, then the B'th cell must have something in it already.

See the FAQ to get a good intuitive feel for what a cell array is.

https://matlab.fandom.com/wiki/FAQ#What_is_a_cell_array.3F

Answer 2

[Image Analyst]

Evidently, it never gets to the line

ZZ{B}(nP) = zAxis{B}(index(1)) ; % We use index(1) because there might be many minimum values 

Set a break point there and see if it does. Otherwise set other breakpoints and follow the execution path of your program by stepping through it:

Debugging in MATLAB | Doug's MATLAB Video Tutorials