Need help in code optimization for distance calculation between points.

Hello everyone. I need to know if there is a way to optimize my code. This part (which concentrates 95% of the computing time) takes about 15 seconds on my computer, and the section that is taking the longest time is calculating the distance (line 20). This time is really large, considering that this is only one region of a larger map and in addition it only equals one iteration of 96. The total time of my program considering all iterations would be 8 hours. Any ideas?

clear;% clc;
load('vars_optimiz.mat','-mat');
time_clc=tic;
polyg='Polyg #1';
fprintf('\n%s (%i:00)',polyg,Dates.H);
n_cell_activ=0;
    alfn(size(DATA_Value,1))=0.0;
    gaman(size(DATA_Value,1))=0.0;
    val_es(size(DATA_Value,1))=0.0;
    sum_m1=0;
    sum_m2=0;
    for x=1:size(mapa,2)
        for y=1:size(mapa)
            sum2_th=0;
            sum2_id=0;
            if mapa(y,x) ~= 0
                n_cell_activ=n_cell_activ+1;
                xx=base_x+x; yy=base_y+y;
                alft=0;
                n_dist=(((DATA_x_s-xx) .^ 2) + ((DATA_y_s-yy) .^ 2)) .^ 0.5; %Distance Calculation
                [val_min_dist,cell_min_dist]=min(n_dist);
                n_dist(n_dist==0)=0.0000000000000001;
                sum_m1=sum_m1+DATA_Value(cell_min_dist);
                alfn=1./(n_dist.^2);
                alft=sum(alfn);
                gaman=alfn/alft;
                val_es=gaman.*DATA_Value;
                val_tot=sum(val_es);
                sum_m2=sum_m2+val_tot;
            end
        end
    end
    med1=sum_m1/n_cell_activ;
    med2=sum_m2/n_cell_activ;
    fprintf(' "%s",%i-%i-%i %i:00 - Med1: %.8f Med2: %.8f Segs: %.4f\n',polyg,Dates.Y,Dates.M,Dates.D,Dates.H,med1,med2,toc(time_clc));

Answer 1

[Jan]

A first step is to clean the loop:

n_dist=(((DATA_x_s-xx) .^ 2) + ((DATA_y_s-yy) .^ 2)); %Distance Calculation
[val_min_dist,cell_min_dist] = min(n_dist);
val_min_dist = sqrt(val_min_dist);
n_dist(n_dist==0) = 1.0e-32;

Replacing ().^0.5 by sqrt() reduces the total runtime by 10% already. The minimum of the squareroot of the values equals teh squareroot of the minium. Because later on you square n_dist again, there is no reason to calculate the squareroot of all values.

This reduces the runtime from 22 to 16 seconds.

Now move the repeated calculations out of the loop. E.g. sum(X / n) for a scalar n and the vector X equals sum(X) / n, but this requires less divisions.

Poofing variables in the workspace by a load without catching the output can impede the JIT acceleration. Therefore the variables are imported explicitely now.

function optimiz
data = load('vars_optimiz.mat','-mat');
base_x     = data.base_x;
mapa       = data.mapa;
base_y     = data.base_y;
Dates      = data.Dates;
DATA_x_s   = data.DATA_x_s;
DATA_y_s   = data.DATA_y_s;
DATA_Value = data.DATA_Value;

time_clc = tic;
polyg = 'Polyg #1';
fprintf('\n%s (%i:00)',polyg, Dates.H);
n_cell_activ = 0;
sum_m1 = 0;
sum_m2 = 0;
for x = 1:size(mapa,2)
    xx = base_x + x;
    cx = (DATA_x_s - xx) .^ 2;
    for y = 1:size(mapa, 1)  % Not 1:size(mapa)
        if mapa(y,x) ~= 0
            n_cell_activ = n_cell_activ+1;
            yy     = base_y + y;
            n_dist = cx + (DATA_y_s - yy) .^ 2;  % Distance Calculation
            [val_min_dist,cell_min_dist] = min(n_dist);
            n_dist(n_dist==0) = 1.0e-32;
            sum_m1  = sum_m1 + DATA_Value(cell_min_dist);
            alfn    = 1 ./ n_dist;
            alft    = sum(alfn);
            val_tot = sum(alfn .* DATA_Value) / alft;
            sum_m2  = sum_m2 + val_tot;
        end
    end
end
med1 = sum_m1 / n_cell_activ;
med2 = sum_m2 / n_cell_activ;
fprintf(' "%s",%i-%i-%i %i:00 - Med1: %.8f Med2: %.8f Segs: %.4f\n', ...
        polyg,Dates.Y,Dates.M,Dates.D,Dates.H,med1,med2,toc(time_clc));

Thsi reduces the run time from 22 to 8.6 seconds. Perhaps this can save further 0.7 seconds:

n_dist(n_dist==0) = 1.0e-32;

can be replaced by:

n_dist = max(n_dist, 1.0e-32);

but this changes the result slightly, when there are some tiny values in n_dist.

Now the loops are cleaner, but the code is not vectorized. Calculating (DATA_y_s - yy) .^ 2 as a matrix might waste time, when mapa(y,x) is zero. Therefore I'm not convinced if a vectorization is useful here. But 3 times faster is a nice start.