Trying to Speed Up Code - Matrix Indexing Help?

Share icon

The background to this problem is that Igor defines atan as between -π/2 and +π/2, and I need the full 360°. So, I wrote a function ("CorrectATan") which not only flips the results when needed to get 0-2π, but also does a little perturbation for those special cases when the value that goes into atan from my work will yield a NaN result (and subsequently mess everything up).

Here's the code ... don't pay attention to all the variables, that's really unimportant for my question:

ThreadSafe Function CorrectATan(x_center_mass_degrees,y_center_mass_degrees,d_planet_radius,rim_lon_temp,rim_lat_temp,rim_lat,rim_lon)
Variable x_center_mass_degrees,y_center_mass_degrees,d_planet_radius
Wave rim_lon_temp,rim_lat_temp,rim_lat,rim_lon
Variable i, d
do
    if(rim_lon[i] < x_center_mass_degrees)
        rim_lon_temp[i] *= -1
        rim_lat_temp[i] *= -1
    else
        if(rim_lon[i] == x_center_mass_degrees) //special case that results in NaN
            d = (sin((rim_lat[i]-y_center_mass_degrees)*pi/180/2))^2 + cos(rim_lat[i]*pi/180)*cos(y_center_mass_degrees*pi/180)*(sin((rim_lon[i]-x_center_mass_degrees)*pi/180/2))^2
            d = 2*atan2(sqrt(d),sqrt(1-d)) * d_planet_radius
            rim_lon_temp[i] = 0
            rim_lat_temp[i] = d * (rim_lat[i]-y_center_mass_degrees)/(((rim_lon[i]+1e-5)-x_center_mass_degrees)*sqrt(1.+(rim_lat[i]-y_center_mass_degrees)^2/((rim_lon[i]+1e-5)-x_center_mass_degrees)^2))
        endif
    endif
    i += 1
while(i < numpnts(rim_lon))
End


This code is fairly slow because it has a do-while loop that goes over every single point, when typically only half the points actually need to be switched over. I thought of doing a Sort and then BinarySearch, but sorts are much slower (O(N*log(N))) than simple iterative searching (O(N)).

After recently learning of the Extract, I thought I'd try that out, if nothing else for the first if-else statement to pull that outside the loop:

Extract/O/INDX rim_lon,index_Temp, (rim_lon < x_center_mass_degrees)
if(numpnts(index_Temp) > 0)
    rim_lon_temp[index_Temp[]] *= -1
    rim_lat_temp[index_Temp[]] *= -1
endif


Problem is, apparently I can't index that way on the left side of the *= sign, and index_Temp may be equal, to, say, {1, 3, 4, 5, 15, 16, 17}. So, I had to do this:

Extract/O/INDX rim_lon,index_Temp, (rim_lon < x_center_mass_degrees)
if(numpnts(index_Temp) > 0)
    do
        rim_lon_temp[index_Temp[i]] *= -1
        rim_lat_temp[index_Temp[i]] *= -1
        i += 1
    while(i < numpnts(index_Temp))
endif


And now I'm at a point where my code takes longer to run than it did before, probably because it still has to do a loop over ~half the points and also do the Extract now.

So ... any ideas? Is it at all possible to index "rim_lon_temp" and "rim_lat_temp" the way I want to? Or is there a better solution?

To give you an idea of time-savings, on code that took 2 minutes to run initially, if I were to remove the CorrectATan part, it runs in 80 seconds. And that's on a small part of my database. So, this is something that could save hours in the long run.
Quote:

atan2(y1, x1 )
The atan2 function returns the angle in radians whose tangent is y1 / x1. Results are in the range -pi to pi.

--Jim Prouty
Software Engineer, WaveMetrics, Inc.

September 7, 2014 at 12:01 am - Permalink

The issue's not quite that simple, and perhaps my function is poorly named (and I just spent the last 15 minutes reading my comment blocks to figure out what I was doing). The problem is: I have points of a circle (lat/lon) in decimal degrees (rim_lat, rim_lon). I know the center (x_center_mass, y_center_mass). I want to convert those rim points in decimal degrees to kilometers given a planet's radius (d_planet_radius, d=double ... carry-over from C).

I'm doing that by using the simpler Great Circle formula from here.

dist = (sin((rim_lat-y_center_mass_degrees)*pi/180/2))^2 + cos(rim_lat*pi/180)*cos(y_center_mass_degrees*pi/180)*(sin((rim_lon-x_center_mass_degrees)*pi/180/2))^2
dist = 2*atan2(sqrt(dist),sqrt(1-dist)) * d_planet_radius


So, I have my distance. Now I need to know the actual angle between the center and my rim points so that I can put the rim points at not only the correct distance from the center, but also the correct angle.

Doing that uses a lot of PITA math (and requires that special case I mentioned in my initial post because you can sometimes get a division by 0):

rim_lon_temp = dist * 1./sqrt(1.+(rim_lat-y_center_mass_degrees)^2/(rim_lon-x_center_mass_degrees)^2)
rim_lat_temp = dist * (rim_lat-y_center_mass_degrees)/((rim_lon-x_center_mass_degrees)*sqrt(1.+(rim_lat-y_center_mass_degrees)^2/(rim_lon-x_center_mass_degrees)^2))


When that's done, I have my angle, but only between ±π/2. So it's not actually an ATan correction (I may have thought it was when I initially programmed it 13 months ago), but rather a solving-for-angles-maybe-±sqrt-issue.

I don't want the overhead of solving a full Vincenty set of equations because that level of accuracy isn't needed, and it would take significantly longer to iterate until a solution converged.

Anyway, that's why I said to ignore the purpose of the algorithm and focus more on the general idea of a way to make that kind of indexing faster, if possible.

September 7, 2014 at 11:38 am - Permalink

Why a do-while with if else if ... end if over all points instead of a direct pass over the wave with implicit test constructions?

Perhaps push the special case function out (though I could not speak to its ThreadSafe compatibility) ...

rim_lon_temp = rim_lon[p] < x_center_mass_degrees ? -rim_lon_temp : rim_lon_temp
rim_lon_temp = rim_lon[p] == x_center_mass_degrees ? rim_lon_temp : 0
rim_lat_temp = rim_lon[p] < x_center_mass_degrees ? -rim_lat_temp : rim_lat_temp
rim_lon_temp = rim_lon[p] == x_center_mass_degrees ? rim_lon_temp : MySpecialATANFunction(rim_lat[p], y_center_mass_degrees, d)

--
J. J. Weimer
Chemistry / Chemical & Materials Engineering, UAHuntsville

September 7, 2014 at 06:58 pm - Permalink

Ahh ... if I had known about such a wonderful iif / ?: mechanism, I might've tried that earlier. My only formal programming class was in high school in 2000-2001.

I also just realized that my special case for the latitude reduces to just d. Very long and complicated way to get all those numbers to equal 1.

September 7, 2014 at 09:45 pm - Permalink