Trying to Speed Up Code - Matrix Indexing Help?
astrostu
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
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
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
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.
--Jim Prouty
Software Engineer, WaveMetrics, Inc.
September 7, 2014 at 12:01 am - Permalink
I'm doing that by using the simpler Great Circle formula from here.
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_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
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 : 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
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