Vectorize my code in azimuth calculation.

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

[Adam]

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]

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

Answer 1

[Guillaume]

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.

[Mani Ahmadian]

Dear Guillaume;

Your answer shocked me, It's great and is very useful for me. I run it and elapsed time is about 0.222 seconds. This time is acceptable for my project and I used your offer in my code.

Another stage is to find all azimuths are between 180+-30 (150<=azimuths and azimuths<=210). I use my code as bellow to find all azimuths are true, it's so fast as I need:

TrueAzimuth=zeros(size(AzimuthPoints));
TrueAzimuth(AzimuthPoints()<=UserAz+AzTol & AzimuthPoints()>=UserAz-AzTol)=1;

As final part I should pick all lat/long values that it's corresponding value in TrueAzimuth is 1. For example if TrueAzimuth(1,3)=1, I pick lat/long values of point#1 and point#3 for future calculations. To select points, I use my simple code as bellow, But it's so slow.

for ii=1:length(lat)
    for jj=1:length(long)
        if TrueAzimuth(ii,jj)==1
            TrueLat=vertcat(TrueLat, lat(ii), lat(jj));
            TrueLong=vertcat(TrueLong, long(ii), long(jj));
        end
    end
end

Please help me to vectorize it, too.

Thanks a lot.

Mani

[Guillaume]

First, TrueAzimuth is just

TrueAzimuth = AzimuthPoints<=UserAz+AzTol & AzimuthPoints>=UserAz-AzTol; %No need to predeclare.

Then, I believe the following will give you exactly the same result as your loop:

TrueLat = [lat1(TrueAzimuth) lat2(TruAzimuth)]';
TrueLat = TrueLat(:);
%same with TrueLong

But the following seems to make more sense to me, if you're going to make further processing:

TrueLat = [lat1(TrueAzimuth) lat2(TrueAzimuth)]; %keep as 2 column matrix

[Mani Ahmadian]

I use your code, but because we're useing TtueAzimuth(consist of 0,1) to index lat1,long1,... I receive an error message: Subscript indices must either be real positive integers or logicals.

How can I solve it?

Thanks so much

Mani

[Guillaume]

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

[Mani Ahmadian]

Thanks, It works correctly.

Answer 2

[Roger Stafford]

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