unrecognized function or variable in parfor

8 visualizaciones (últimos 30 días)
Maurilio Matracia
Maurilio Matracia el 17 de Sept. de 2021
Comentada: Image Analyst el 20 de Sept. de 2021
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
  2 comentarios
Geoff Hayes
Geoff Hayes el 17 de Sept. de 2021
Maurilio - where is ZZ defined?
Maurilio Matracia
Maurilio Matracia el 17 de Sept. de 2021
Inside Montecarlo integration section, almost at the end

Iniciar sesión para comentar.

Respuestas (2)

Raymond Norris
Raymond Norris el 17 de Sept. de 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 comentarios
Maurilio Matracia
Maurilio Matracia el 19 de Sept. de 2021
Thank you very much. Actually I just had to initialize ZZ and zAxis to make it work, do you know why?
Image Analyst
Image Analyst el 20 de Sept. de 2021
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.

Iniciar sesión para comentar.


Image Analyst
Image Analyst el 17 de Sept. de 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:

Categorías

Más información sobre Call Python from MATLAB en Help Center y File Exchange.

Community Treasure Hunt

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

Start Hunting!

Translated by