[parfor] Error: The variable ViewF in a parfor cannot be classified.

Bonjour à tous,

J'ai un petit souci avec l'utilisation de parfor.

Quand je fais tourner le code ci-dessous je n'ai aucun problème :

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
 
tic
time_zone = 4 ;                                                            %Time zone (i.e. GMT+) ;
dst = 0 ;                                                                     %Daylight saving time ;                     
 
Lmax = 55.8272 ;                                                           %Goegraphic area for calculation ;
Lmin = 55.2168 ; 
lmax = -20.8722 ; 
lmin = -21.3891 ; 
 
step = 0.1 ; 
LL = Lmin:step:Lmax ; 
ll = lmin:step:lmax ;
 
R_az = 45  ;                                                               %Azimuth angle resolution [°] ;
R_di = 0.005 ;                                                              %Resolution a direction defined by any azimuth angle [°] ;
 
R_az = R_az * pi / 180 ;                                                   %Conversion degre to radian ; 
 
matlabpool('open',2);
 
parfor n = 1:length(LL)
    for m = 1:length(ll)
 
        lat = ll(m) ;                                                      %Geopoint M(lat,long) ;
        lon = LL(n) ;
 
        display(lat);
        display(lon); 
        display( (now - time1)*86400 ) ;  
 
        altM = interp2( lat_long_alt_Reu_STRM(:,:,2),...                   %Altitude of M [m];
                               lat_long_alt_Reu_STRM(:,:,1),...
                               lat_long_alt_Reu_STRM(:,:,3),...
                               lon , lat ,'cubic') ;
 
        if isnan( altM ) 
 
            altM = 0 ;             
 
        end                       
 
        azi = pi:-R_az:-pi ;
        Hz = zeros(1,length(azi));
 
        for i = 1:length(azi)                                              %Azimuth step angle vector ;  
            k = 1 ; 
            r = 0 ; 
 
            X = lon ; 
            Y = lat ; 
 
            while  lmin <= Y && Y <= lmax && Lmin <= X && X <= Lmax ;
 
                r = r + R_di ;                                             %Geopoint N(X,Y) ;
                X = lon + r * cos(azi(i)+pi/2) ; 
                Y = lat + r * sin(azi(i)+pi/2) ;
 
                altN = interp2( lat_long_alt_Reu_STRM(:,:,2),...           %Altitude of N [m] ;
                               lat_long_alt_Reu_STRM(:,:,1),...
                               lat_long_alt_Reu_STRM(:,:,3),...
                               X , Y ,'cubic') ; 
 
                if isnan( altN ) 
 
                    altN = 0 ; 
 
                end    
 
 
                altMN = altN - altM ;                                      %Altitude difference between M and N [m] ;
 
                distMN = vdist( lat , lon , Y , X ) ;                      %Distance between M and N [m] ;
 
                Ele(n,k) = atan( altMN / distMN ) * 180 / pi ;               %Elevation of N from M ;
                k = k + 1 ;
 
            end   
 
            Hz(i) = max(Ele(n,:)) ;                                             %Overshagin of M ;
            Ele = zeros(1) ;
        end
 
        azi(1) = [] ; 
        Hz(1) = [] ;
 
        %Sky View Factor calculation: 
        TotalSky = 0 ;
        SkyView = 0 ;
        for i = 1:length(Hz) 
 
            for j = 1:0.1:90 ; 
 
                TotalSky = TotalSky + 1 ;   
 
                if j > Hz(i) 
 
                    SkyView = SkyView + 1 ;  
 
                end
 
            end
 
        end 
 
        SVF = SkyView / TotalSky ; 
 
        %Sun Beam (and weighted) View Factor calculation: 
        TotalSunPath = 0 ; 
        SunPathView = 0 ;  
        TotalBeamEnergy = 0 ; 
        BeamEnergyView = 0 ; 
 
        for j=1:365
            for i=1:1440  
 
                %Note: i represents the minute of the day and j represents the
                %day of the year. 
 
                %The equation of time (E):
                B = ( j - 1 ) * 2 * pi / 365;
                E = 229.18 * ( 7.5e-5 + 1.868e-3 * cos(B) - 3.2077e-2 * sin(B)- 1.4615e-2 * cos(2*B) - 4.089e-2 * sin(2*B) );
                %Convert to west longitude:
                if lon>0
                    longitude_standard=360-(time_zone+dst)*15;
                    longitude_in_west_degree=360-lon; 
                else
                    longitude_standard= -(time_zone+dst)*15;
                    longitude_in_west_degree= -lon;
                end
 
                %Difference in minutes between solar time and local time
                solartime_minus_localtime = 4 *( longitude_standard - longitude_in_west_degree )+ E;
 
                %Solar time (minute):
                solar_time= i + solartime_minus_localtime;
                if solar_time > 1440 
                solar_time = solar_time - 1440 ; 
                end    
 
                if solar_time <= 0 
                    solar_time = 1440 - solar_time ; 
                end 
 
                %Decimal solar hour
                solar_hour = solar_time/60 ; 
 
                %Solar angle calculation:
                %Hour angle (omega):
                omega = (solar_hour - 12) * 15 * pi / 180 ;
 
                %Declination angle (delta):
                delta = 23.45 * pi/180 * sin (2*pi*(j+284)/365); 
 
                %Solar elevation angle (alpha_s):
                alpha_s = asin( cos(omega) * cos(delta) * cos(lat*pi/180) + sin(delta) * sin(lat*pi/180) );
                alpha_s_degree = 180 * alpha_s/pi;
 
                % Solar azimuth angle (gamma_s):
                if lat >= 0
                    gamma_s = atan2( sin(omega) , cos(omega)*sin(lat*pi/180)-tan(delta)*cos(lat*pi/180) ) + pi ;
                else
                    gamma_s = atan2( sin(omega) , cos(omega)*sin(lat*pi/180)-tan(delta)*cos(lat*pi/180) ) ;
 
                    if gamma_s > 0
                        gamma_s = gamma_s - pi ;
                    else
                        gamma_s = gamma_s + pi ; 
                    end
 
                end
 
                gamma_s_degree = gamma_s * 180 / pi ;
 
                if alpha_s_degree > 90 
 
                    theta_s = ( alpha_s_degree - 90 ) * pi / 180 ;
 
                else
 
                    theta_s =  ( 90 - alpha_s_degree ) * pi / 180 ; 
 
                end 
 
                if alpha_s > 0 
 
                    RavR2 = 1.0001 + 0.034221 * cos(B) + 0.001280 * sin(B) + 0.000719 * cos(2*B) + 0.000077 * sin(2*B); 
                    %Note: RavR2 is the ratio beteween the mean sun-earth distance 
                    %(Rav) and the the actual sun-earth distance depending on the day of the year 
 
                    I0 = 1367 * RavR2 ; %Extraterrestrial radiation in W.m-2
 
                    AM = 1 / cos(theta_s) ;  
 
                    BNI = I0*0.7^(AM^0.678) ; 
 
                    TotalSunPath = TotalSunPath + 1 ; 
                    TotalBeamEnergy = TotalBeamEnergy + BNI ; 
 
                    ind = find( azi <= gamma_s ) ; 
 
                    if alpha_s_degree > Hz(ind(1))
 
                        SunPathView = SunPathView + 1 ; 
                        BeamEnergyView = BeamEnergyView + BNI ; 
 
                    end
 
                end
 
            end
        end
 
        SBVF = SunPathView / TotalSunPath ; 
        SBWVF = BeamEnergyView / TotalBeamEnergy ; 
 
 
    end
end
 
matlabpool('close');
 
toc

Le calcul s'effectue bien sur deux cœurs et à chaque fois que je passe dans les boucle n et m je calcul bien les variables : lat , lon , SVF , SBVF et SBWF.
Maintenant je souhaite récupérer ces variables pour chaque couple (n,m) calculé. Pour cela j'introduis le tableau ViewF, mais Matlab me renvoi l'erreur suivante :

Error: The variable ViewF in a parfor cannot be classified.

Voici plus bas comment j'utilise ViewF:

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
%*************************************************************************%
tic
time1 = now ; 
disp('processing...')
 
time_zone = 4 ;                                                            %Time zone (i.e. GMT+) ;
dst = 0 ;                                                                  %Daylight saving time ;                     
 
Lmax = 55.8272 ;                                                           %Goegraphic area for calculation ;
Lmin = 55.2168 ; 
lmax = -20.8722 ; 
lmin = -21.3891 ; 
 
step = 0.1 ; 
LL = Lmin:step:Lmax ; 
ll = lmin:step:lmax ;
 
R_az = 45  ;                                                               %Azimuth angle resolution [°] ;
R_di = 0.005 ;                                                              %Resolution a direction defined by any azimuth angle [°] ;
 
R_az = R_az * pi / 180 ;                                                   %Conversion degre to radian ; 
 
 
spmd
    ViewF = zeros(length(LL), length(ll) ,5 ); 
end
 
 
matlabpool('open',2);
 
parfor n = 1:length(LL)
    for m = 1:length(ll)
 
        lat = ll(m) ;                                                      %Geopoint M(lat,long) ;
        lon = LL(n) ;
 
        display(lat);
        display(lon); 
        display( (now - time1)*86400 ) ;  
 
        altM = interp2( lat_long_alt_Reu_STRM(:,:,2),...                   %Altitude of M [m];
                               lat_long_alt_Reu_STRM(:,:,1),...
                               lat_long_alt_Reu_STRM(:,:,3),...
                               lon , lat ,'cubic') ;
 
        if isnan( altM ) 
 
            altM = 0 ;             
 
        end                       
 
        azi = pi:-R_az:-pi ;
        Hz = zeros(1,length(azi));
 
        for i = 1:length(azi)                                              %Azimuth step angle vector ;  
            k = 1 ; 
            r = 0 ; 
 
            X = lon ; 
            Y = lat ; 
 
            while  lmin <= Y && Y <= lmax && Lmin <= X && X <= Lmax ;
 
                r = r + R_di ;                                             %Geopoint N(X,Y) ;
                X = lon + r * cos(azi(i)+pi/2) ; 
                Y = lat + r * sin(azi(i)+pi/2) ;
 
                altN = interp2( lat_long_alt_Reu_STRM(:,:,2),...           %Altitude of N [m] ;
                               lat_long_alt_Reu_STRM(:,:,1),...
                               lat_long_alt_Reu_STRM(:,:,3),...
                               X , Y ,'cubic') ; 
 
                if isnan( altN ) 
 
                    altN = 0 ; 
 
                end    
 
 
                altMN = altN - altM ;                                      %Altitude difference between M and N [m] ;
 
                distMN = vdist( lat , lon , Y , X ) ;                      %Distance between M and N [m] ;
 
                Ele(n,k) = atan( altMN / distMN ) * 180 / pi ;               %Elevation of N from M ;
                k = k + 1 ;
 
            end   
 
            Hz(i) = max(Ele(n,:)) ;                                             %Overshagin of M ;
            Ele = zeros(1) ;
        end
 
        azi(1) = [] ; 
        Hz(1) = [] ;
 
        %Sky View Factor calculation: 
        TotalSky = 0 ;
        SkyView = 0 ;
        for i = 1:length(Hz) 
 
            for j = 1:0.1:90 ; 
 
                TotalSky = TotalSky + 1 ;   
 
                if j > Hz(i) 
 
                    SkyView = SkyView + 1 ;  
 
                end
 
            end
 
        end 
 
        SVF = SkyView / TotalSky ; 
 
        %Sun Beam (and weighted) View Factor calculation: 
        TotalSunPath = 0 ; 
        SunPathView = 0 ;  
        TotalBeamEnergy = 0 ; 
        BeamEnergyView = 0 ; 
 
        for j=1:365
            for i=1:1440  
 
                %Note: i represents the minute of the day and j represents the
                %day of the year. 
 
                %The equation of time (E):
                B = ( j - 1 ) * 2 * pi / 365;
                E = 229.18 * ( 7.5e-5 + 1.868e-3 * cos(B) - 3.2077e-2 * sin(B)- 1.4615e-2 * cos(2*B) - 4.089e-2 * sin(2*B) );
                %Convert to west longitude:
                if lon>0
                    longitude_standard=360-(time_zone+dst)*15;
                    longitude_in_west_degree=360-lon; 
                else
                    longitude_standard= -(time_zone+dst)*15;
                    longitude_in_west_degree= -lon;
                end
 
                %Difference in minutes between solar time and local time
                solartime_minus_localtime = 4 *( longitude_standard - longitude_in_west_degree )+ E;
 
                %Solar time (minute):
                solar_time= i + solartime_minus_localtime;
                if solar_time > 1440 
                solar_time = solar_time - 1440 ; 
                end    
 
                if solar_time <= 0 
                    solar_time = 1440 - solar_time ; 
                end 
 
                %Decimal solar hour
                solar_hour = solar_time/60 ; 
 
                %Solar angle calculation:
                %Hour angle (omega):
                omega = (solar_hour - 12) * 15 * pi / 180 ;
 
                %Declination angle (delta):
                delta = 23.45 * pi/180 * sin (2*pi*(j+284)/365); 
 
                %Solar elevation angle (alpha_s):
                alpha_s = asin( cos(omega) * cos(delta) * cos(lat*pi/180) + sin(delta) * sin(lat*pi/180) );
                alpha_s_degree = 180 * alpha_s/pi;
 
                % Solar azimuth angle (gamma_s):
                if lat >= 0
                    gamma_s = atan2( sin(omega) , cos(omega)*sin(lat*pi/180)-tan(delta)*cos(lat*pi/180) ) + pi ;
                else
                    gamma_s = atan2( sin(omega) , cos(omega)*sin(lat*pi/180)-tan(delta)*cos(lat*pi/180) ) ;
 
                    if gamma_s > 0
                        gamma_s = gamma_s - pi ;
                    else
                        gamma_s = gamma_s + pi ; 
                    end
 
                end
 
                gamma_s_degree = gamma_s * 180 / pi ;
 
                if alpha_s_degree > 90 
 
                    theta_s = ( alpha_s_degree - 90 ) * pi / 180 ;
 
                else
 
                    theta_s =  ( 90 - alpha_s_degree ) * pi / 180 ; 
 
                end 
 
                if alpha_s > 0 
 
                    RavR2 = 1.0001 + 0.034221 * cos(B) + 0.001280 * sin(B) + 0.000719 * cos(2*B) + 0.000077 * sin(2*B); 
                    %Note: RavR2 is the ratio beteween the mean sun-earth distance 
                    %(Rav) and the the actual sun-earth distance depending on the day of the year 
 
                    I0 = 1367 * RavR2 ; %Extraterrestrial radiation in W.m-2
 
                    AM = 1 / cos(theta_s) ;  
 
                    BNI = I0*0.7^(AM^0.678) ; 
 
                    TotalSunPath = TotalSunPath + 1 ; 
                    TotalBeamEnergy = TotalBeamEnergy + BNI ; 
 
                    ind = find( azi <= gamma_s ) ; 
 
                    if alpha_s_degree > Hz(ind(1))
 
                        SunPathView = SunPathView + 1 ; 
                        BeamEnergyView = BeamEnergyView + BNI ; 
 
                    end
 
                end
 
            end
        end
 
        SBVF = SunPathView / TotalSunPath ; 
        SBWVF = BeamEnergyView / TotalBeamEnergy ; 
 
        ViewF(n,m,1) = lat ;  
        ViewF(n,m,2) = lon ; 
        ViewF(n,m,3) = SVF ;
        ViewF(n,m,4) = SBVF ; 
        ViewF(n,m,5) = SBWVF ; 
 
    end
end
 
matlabpool('close');
 
toc

Je ne comprends pas ce qui pose problème.

Merci par avance

je te conseille ce lien : http://www.mathworks.fr/fr/help/dist...ed-topics.html

en gros l'idée c'est que dans ton tableau tu peux pas mettre en même temps ta variable de boucle et d'autre variables, je cite :

Citation:

With i as the loop variable, the A variables shown here on the left are not sliced; those on the right are sliced:

Not sliced Sliced
A(i+f(k),j,:,3) A(i+k,j,:,3)
A(i,20:30,end) A(i,:,end)
A(i,:,s.field1) A(i,:,k)

When you use other variables along with the loop variable to index an array, you cannot set these variables inside the loop. In effect, such variables are constant over the execution of the entire parfor statement. You cannot combine the loop variable with itself to form an index expression.

du coup tu peux essayer (je sais pas si ça va marcher, j'ai pas la license pour la parallel computing toolbox) :
-soit d'inverser l'ordre de ton for et de ton parfor
Code:

1 2 3 4 5 6 for m = 1:length(ll) parfor n = 1:length(LL) .... end end
-soit de sauver un "vecteur temporaire" pour le mettre dans ton tableau plus tard
Code:

1 2 3 4 5 6 7 8 9 10 11 parfor n = 1:length(LL) for m = 1:length(ll) ViewF_vect(m,1) = lat ; ViewF_vect(m,2) = lon ; ViewF_vect(m,3) = SVF ; ViewF_vect(m,4) = SBVF ; ViewF_vect(m,5) = SBWVF ; end ViewF(n,:,:) = ViewF_vect; end
je priviligierai personellement la première option, je ne suis pas très sur que la seconde marche.

Tu ne peux pas utiliser une variable Composite issue d'un bloc spmd dans une boucle parfor

Dans ton cas, tu peux simplement initialiser ViewF avant d'entrer dans la boucle parfor :

Code:

ViewF = zeros(length(LL), length(ll) ,5 );

Salut et merci pour vos conseils.

J'ai testé plusieurs choses et voici ce qui a fonctionné :

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
 
ViewF = zeros(length(LL), length(ll) ,5 ); 
 
parfor n = 1:length(LL)
 
    ViewF_vect = zeros( length(ll) ,5 ); 
 
    for m = 1:length(ll)
 
        [...]
 
        ViewF_vect(m,1) = lat ;  
        ViewF_vect(m,2) = lon ; 
        ViewF_vect(m,3) = SVF ;
        ViewF_vect(m,4) = SBVF ; 
        ViewF_vect(m,5) = SBWVF ; 
 
     end 
 
     ViewF(n,:,:) =  ViewF_vect;
 
end

Merci

Ceci devrait également fonctionner :

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
nll = length(ll);
 
ViewF = zeros(length(LL),nll ,5 ); 
 
parfor n = 1:length(LL)
 
    for m = 1:nll
 
        [...]
 
        ViewF(n,m,1) = lat ;  
        ViewF(n,m,2) = lon ; 
        ViewF(n,m,3) = SVF ;
        ViewF(n,m,4) = SBVF ; 
        ViewF(n,m,5) = SBWVF ; 
 
     end 
 
end

Voir "Limitations of Nested for-Loops" sur cette page : Programming Considerations