how can I reduce the computational time of a single for loop in a Matlab code

1 view (last 30 days)
I have a problem with the following for. it is too time consuming. could anyone help me with it?thank in advance for your consideration and time.
Let me give you some more details about the aim of the code. I have observations with the number of 'd' that measured in 'day-of-year'. so each observations occured in a specific 'day-of-year'.
I want to collocate each of these observations (10 million data) with 'VWC', 'SM_R' that measured in coordinates (Lat, Lon) and (LAT_R, LON_R) in days called 'DOY' and 'DOY-R' in code, respectivelly. see the associate part of code:
%%%%%Matlab Code%%%%%%
Lon_d=Lon(DOY==day_of_year(d));
Lat_d=Lat(DOY==day_of_year(d));
VWC_SMAP_d=VWC(DOY==day_of_year(d));
LAT_R_d=LAT_R(DOY_R==day_of_year(d));
LON_R_d= LON_R(DOY_R==day_of_year(d));
SM_R_9_d=SM_R(DOY_R==day_of_year(d));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
to achieve this goal, each observation with coordinate 'lon_cy(d), lat_cy(d)' was collocated with observations limited to 'step (18/110 degree) *step (18/110 degree)' grids with the center of it (lon_cy(d), lat_cy(d)):
%%%%%%%Matlab Code%%%%%%%%%%%%%%%%%%%%
lon_smap_9_index= find(LON_R_d <=lon_cy(d)+step & LON_R_d>=lon_cy(d)-step);
lat_smap_9_index= find(LAT_R_d <=lat_cy(d)+step & LAT_R_d>=lat_cy(d)-step);
intsec_smap_9=intersect(lon_smap_9_index,lat_smap_9_index);
VWC_SMAP_g_9=VWC_R_9_d(intsec_smap_9)
lon_smap_36_index= find(Lat_d <=lat_cy (d)+step & Lat_d>=lat_cy(d)-step);
lat_smap_36_index= find(Lon_d <=lon_cy(d)+step & Lon_d>=lon_cy(d)-step);
intsec_smap_36=intersect(lon_smap_36_index,lat_smap_36_index);
SM_g_36=SM_d(intsec_smap_36);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
two additional ancillary datasets were aimed to colloated with the d-dimension observations for each day-of-year:
here I just need to find the nearest ancillary datasets to each observation. that is why I used 'knnsearch' function.
%%%%%%%%%%%%%%%%%%%Matlab%%%%%%%%%%%%%%
point_1=[lat_cy(d) lon_cy(d)];
ptCloud_1=[latitu longi];
idx_land = knnsearch(ptCloud_1,point_1);
ptCloud_2=[clay(:,1) clay(:,2)];
idx_clay = knnsearch(ptCloud_2,point_1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
and finally I will store each collocated observations in the following way and save it
%%%%%%%%%%%%%%%%%%%%Matlab%%%%%%%%%%%%%%%
FT=[day_of_year(d) lon_cy(d) lat_cy(d) inc_cy(d) ref_cy(d) mean(VWC_SMAP_g_9) mean(roughness_g_9) mean(SM_g_9) mean(VWC_SMAP_g_36) mean(roughness_g_36) mean(TM_d_36) mean(SM_g_36) clay(idx_clay,3) land_cover(idx_land)];
  8 Comments
Walter Roberson
Walter Roberson on 8 Apr 2023
Lon_d=Lon(DOY==day_of_year(d));
Lat_d=Lat(DOY==day_of_year(d));
VWC_SMAP_d=VWC(DOY==day_of_year(d));
roughness_d=roughness(DOY==day_of_year(d));
TM_d=TM(DOY==day_of_year(d));
SM_d=SM_36(DOY==day_of_year(d));
LAT_R_d=LAT_R(DOY_R==day_of_year(d));
LON_R_d= LON_R(DOY_R==day_of_year(d));
VWC_R_9_d=VWC_R(DOY_R==day_of_year(d));
SM_R_9_d=SM_R(DOY_R==day_of_year(d));
ROUGH_R_9_d=ROUGH_R(DOY_R==day_of_year(d));
You can pre-calculate those for each different value 1 to 365 (or 366 if appropriate), and index them by day_of_year(d), instead of calculating them every time.
Mina Mino
Mina Mino on 8 Apr 2023
Edited: Mina Mino on 8 Apr 2023
thanks for your reply. I rewrite the code in the following form, the elapsed time decreased from 2.3 to 0.3 second for each day-of-year. However, it also takes a days and even months for 12 million or 20 milion observations. I think I have to change the way of coding. several days ago one of our friends here recommend a way of thinking for another code which reduce the computational time from months to a few second. Do you have any idea?
%%
FT=zeros(1,14);step=18/110;
SMAP_36=[DOY Lon Lat VWC roughness TM SM_36 ];
SMAP_9=[DOY_R LON_R LAT_R VWC_R ROUGH_R SM_R];
parfor d=8350381:10350380
tic
SMAP_36_d=SMAP_36(SMAP_36(:,1)==day_of_year(d),:);
SMAP_9_d=SMAP_9(SMAP_9(:,1)==day_of_year(d),:);
if (~isempty(SMAP_36_d) && ~isempty(SMAP_9_d))
lon_smap_9_index= find(SMAP_36_d(:,2) <=lon_cy(d)+step & SMAP_36_d(:,2)>=lon_cy(d)-step);
lat_smap_9_index= find(SMAP_36_d(:,3) <=lat_cy(d)+step & SMAP_36_d(:,3)>=lat_cy(d)-step);
intsec_smap_9=intersect(lon_smap_9_index,lat_smap_9_index);
lon_smap_36_index= find(SMAP_9_d(:,3) <=lat_cy (d)+step & SMAP_9_d(:,3)>=lat_cy(d)-step);
lat_smap_36_index= find(SMAP_9_d(:,2) <=lon_cy(d)+step & SMAP_9_d(:,2)>=lon_cy(d)-step);
intsec_smap_36=intersect(lon_smap_36_index,lat_smap_36_index);
if (~isempty(intsec_smap_36)&& ~isempty(intsec_smap_9))
VWC_SMAP_g_9=SMAP_9(intsec_smap_9,4);roughness_g_9=SMAP_9(intsec_smap_9,5);SM_g_9=SMAP_9(intsec_smap_9,6);
VWC_SMAP_g_9(isnan(VWC_SMAP_g_9))=[];
roughness_g_9(isnan(roughness_g_9))=[];
SM_g_9(isnan(SM_g_9))=[];
VWC_SMAP_g_36=SMAP_36(intsec_smap_36,4);roughness_g_36=SMAP_36(intsec_smap_36,5);TM_d_36=SMAP_36(intsec_smap_36,6) ;SM_g_36=SMAP_36(intsec_smap_36,7);
VWC_SMAP_g_36(isnan(VWC_SMAP_g_36))=[];
roughness_g_36(isnan(roughness_g_36))=[];
TM_d_36(isnan(TM_d_36))=[];
SM_g_36(isnan(SM_g_36))=[];
point_1=[lat_cy(d) lon_cy(d)];
ptCloud_1=[latitu longi];
idx_land = knnsearch(ptCloud_1,point_1);
ptCloud_2=[clay(:,1) clay(:,2)];
idx_clay = knnsearch(ptCloud_2,point_1);
FT=[day_of_year(d) lon_cy(d) lat_cy(d) inc_cy(d) ref_cy(d) mean(VWC_SMAP_g_9) mean(roughness_g_9) mean(SM_g_9) mean(VWC_SMAP_g_36) mean(roughness_g_36) mean(TM_d_36) mean(SM_g_36) clay(idx_clay,3) land_cover(idx_land)];
S=FT;
% parsave(sprintf('output%d.mat', d),S );
end
end
toc
end

Sign in to comment.

Answers (1)

埃博拉酱
埃博拉酱 on 9 Apr 2023
Edited: 埃博拉酱 on 9 Apr 2023
It seems that you are calling knnsearch twice in each loop to search for a single point in the point cloud. However, knnsearch can actually search for multiple points at once. It is recommended that you try to move knnsearch out of the loop and use the vectorized syntax of knnsearch to calculate all idx_land and idx_clay at once.
Based on my experience, I feel that the bottleneck of your code's speed is this knnsearch. If it's not, please reply to me again.

Categories

Find more on Numeric Types in Help Center and File Exchange

Products


Release

R2021a

Community Treasure Hunt

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

Start Hunting!