Vectorize my code in azimuth calculation.

1 view (last 30 days)
Hi I have lots of points(in Latitude & Longitude) and I want to calculate the azimuth of each pairs of them. I use code as attached file 'AzimuthCalc.m'. when I run my function, it's so time consuming and it's elapsed time is 109.371341 seconds. I'm sure the problem caused by my non-vectorize code.
So, would you please check my function and help me to improve it?
Thanks a lot Mani
  2 Comments
Adam
Adam on 23 Sep 2014
Run:
profile on
your code
profile off
profile viewer
and find out where your code is slow. No point in guessing when you can use facts!
Mani Ahmadian
Mani Ahmadian on 23 Sep 2014
Dear Adam,
Thanks so much. I guess the problem is here(nested for):
for ii=1:length(lat)
for jj=ii+1:length(lat)
AzimuthPoints(ii,jj)=azimuth('rh',lat(ii,1),long(ii,1),lat(jj,1),long(jj,1));
end
end
If I can change my simple code to vectorize one, problem will solve. So, I'm trying to find a way to do that, so plz help me.
Thanks
Mnai

Sign in to comment.

Accepted Answer

Guillaume
Guillaume on 23 Sep 2014
I don't have the mapping toolbox, so I'm just relying on the documentation to write my answer.
For a start, you can replace the inner loop with:
AzimuthPoints(ii,ii+1:end)=azimuth('rh',lat(ii),long(ii),lat(ii+1:end),long(ii+1:end));
I believe the following may work to replace both loops (and the AzimuthPoints declaration)
[lat1, lat2] = meshgrid(lat);
[long1, long2] = meshgrid(long);
AzimuthPoints = azimuth('rh', lat1, long1, lat2, long2);
The downside of this method is that it's going to calculate the whole AzimuthPoints matrix, not just the upper triangle, so it may actually be slower.
You could do:
lat1 = triu(lat1, 1);
%same with lat2, long1 and long2
to replace the lower triangle and main diagonal with zero, on the assumption that calculating
azimuth('rh', 0, 0, 0, 0)
is faster.
  5 Comments
Guillaume
Guillaume on 24 Sep 2014
TrueAzimuth should be 0s and 1s but of type logical. If you used my code to generate it, it shoul already be logical, but if not:
TrueAzimuth = logical(TrueAzimuth);

Sign in to comment.

More Answers (1)

Roger Stafford
Roger Stafford on 24 Sep 2014
You could try this. It calls on 'azimuth' in the same sequence as in your for-loop code, but it is vectorized. Whether it is faster depends on whether the 'tril' and 'repmat' calls take less time than the overhead in the for-loops and individual 'azimuth' calls in your code.
n = length(lat);
p = tril(repmat(1:n,n,1),-1);
p = p(p>0);
q = tril(repmat((1:n)',1,n),-1);
q = q(q>0);
AzimuthPoints = zeros(n);
AzimuthPoints(p+n*(q-1)) = azimuth('rh',lat(p),long(p),lat(q),long(q));

Community Treasure Hunt

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

Start Hunting!