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

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)];

[Dyuman Joshi]

Using logical indexing instead of find() would be a start.

[Mina Mino]

thanks for your reply. I changed it in the following form. But is does not improve significantly. could you please give more tips? Thanks in advance.

%%
FT=zeros(1,14);step=18/110;
parfor  d=8350381:10350380 
    
    
    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));
    
    if (~isempty(ROUGH_R_9_d) && ~isempty(SM_d))
        
        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);
        
        
        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); 
        
        
        VWC_SMAP_g_9=VWC_R_9_d(intsec_smap_9);roughness_g_9=ROUGH_R_9_d(intsec_smap_9);SM_g_9=SM_R_9_d(intsec_smap_9);
        
        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=VWC_SMAP_d(intsec_smap_36);roughness_g_36=roughness_d(intsec_smap_36);TM_d_36=TM_d(intsec_smap_36) ;SM_g_36=SM_d(intsec_smap_36);
        
        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);
        
        
        if (~isempty(intsec_smap_36)&& ~isempty(intsec_smap_9)) 
            
            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
end

[Walter Roberson]

day_of_year suggests something that would have at most 366 different values. A lot of what is calculated is the same for each d value that has the same day_of_year. That suggests you should precalculate those variables for each distinct day_of_year

[Mina Mino]

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)];

[Walter Roberson]

How many unique day_of_year values are there?

[Mina Mino]

all day-of-year,DOY, and DOY-R range from 1 to 365.

[Walter Roberson]

    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]

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

Answer 1

[埃博拉酱]

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.