MATLAB Answers

unrecognized function or variable in parfor

16 views (last 30 days)
Maurilio Matracia
Maurilio Matracia on 17 Sep 2021
Commented: Image Analyst on 20 Sep 2021 at 2:04
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

Answers (2)

Raymond Norris
Raymond Norris on 17 Sep 2021
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.
  2 Comments
Image Analyst
Image Analyst on 20 Sep 2021 at 2:04
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.

Sign in to comment.


Image Analyst
Image Analyst on 17 Sep 2021
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:

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by