How to increase the performance of a cycle for in its use in large matrices?

Dear all,

I have four matrices NewZ(1x20), Longitude(1000x1), Latitude(1000x1) and depthgrid(1000x1000). Note that depthgrid represents the depths in the spatial grid (longitude, latitude). Note also that B is a 1000 x 1000 x 20 matrix that corresponds to the SPL values in the spatial grid along 20 different depths.

I need to check every point of the depthgrid in order to calculate the SPL mean of the points along the longitude x latitude points that are lower than a specific depth (in NewZ matrix).

This is a code I wrote for this purpose:

A=zeros(size(depthgrid));

for i=1:length(longitude);

for j=1:length(latitude);

FINDLOWERS=find(-NewZ<depthgrid(i,j));

BMEDIA=mean(B(:,:,min(FINDLOWERS):max(FINDLOWERS)),3); %To calculate the mean of a cat matrix

A(i,j)=(BMEDIA(i,j)); %place the point into the matrix A

end

end

This works, however it takes me six hours to do the calculations. So, my question is: is it possible to perform a faster checking? Two consecutive for cycles are very uneficient. What should I do?

Thank you very much for your help.

Answer 1

[Jan]

Why do you calculate the mean of the complete BMEDIA matrix, if you use the output for the element (i,j) only?

A = zeros(size(depthgrid));
for i = 1:length(longitude)
    for j = 1:length(latitude)
        F      = find(-NewZ<depthgrid(i,j));      
        A(i,j) = mean(B(i,j, min(F):max(F)),3);
    end
end

mean() is not really efficient. Try this:

A = zeros(size(depthgrid));
for i = 1:length(longitude)
    for j = 1:length(latitude)
        F      = find(-NewZ<depthgrid(i,j));
        F1     = F(1);
        F2     = F(numel(F));
        A(i,j) = sum(B(i,j, F1:F2), 3) / (F2 - F1 + 1);
    end
end

Alternatively:

A  = zeros(size(depthgrid));
N  = zeros(size(depthgrid));
BB = reshape(BB, [], size(B, 3));
nNewZ = -NewZ;
for k = 1:numel(NewZ) - 1
   index    = (nNewZ(k) <= depthgrid < nNewZ(k+1));
   A(index) = A(index) + BB(index, k);
   N(index) = N(index) + 1;
end
A = A ./ N;

There is an indexing problem in this version, therefore I assume the result is not correct yet. Without matching input arguments I cannot debug it currently. Try to fix this code by your own.

[Ricardo Duarte]

What a sorcery.... This is brilliant... Thank you very much! This solved completely the problem.

Another question,

What should I do, if instead the mean I would like to calculate the percentile, 50?

Again, thank you very much Jan.

[Ricardo Duarte]

Done! Thanks again!