function or math problem

[Dana]

I'm working with two scripts, trying to put the functionality of script 2 into script 1. Script 2 reads a .gpx (GPS track) file containing track points with latitude and longitude information, looks for stops, and compares the stop points with another database file containing known locations, finding the closest match. Both scripts contain multiple UDF's. When I wrote script 2 with this in mind, I wrote script 2 with nothing but UDFs; the first thing script 2 does (after declaring some global variables) is to call the first UDF, which is the main routine, which itself calls other UDFs.

I then cut and pasted the entire contents of script 2 into script 1, checking that there are no variable name conflicts. The only changes I made to the script 2 functions were to return back to the main script 1 GUI rather than exiting at the end Depending on user interaction, the UDF which was the main part of script 2 is called, which then calls the other UDFs that were in script 2 and are now in script 1. The problem is, it doesn't work.

I don't expect people to wade through all of both scripts, but the problem seems to originate in one of the smaller UDFs from script 2. This function takes two pairs of latitude and longitude, and calculates the distance between the two:

Func _LLDist($lat1, $lon1, $lat2, $lon2)
;using Sinnott's formula for great circle distance between two points
;first convert degrees to radians for trig functions, using ACos(-1) to derive pi
ConsoleWrite("1197 " & $lat1 & ", " & $lon1 & "; " & $lat2 & ", " & $lon2)
$lat1 = $lat1 * (ACos(-1) / 180)
$lon1 = $lon1 * (ACos(-1) / 180)
$lat2 = $lat2 * (ACos(-1) / 180)
$lon2 = $lon2 * (ACos(-1) / 180)
Local $R = 3440.064794816415 ; Earth radius in nautical miles (6371000 meters)
Local $dlon = $lon2 - $lon1
Local $dlat = $lat2 - $lat1
Local $a = (Sin($dlat / 2)) ^ 2 + Cos($lat1) * Cos($lat2) * (Sin($dlon / 2)) ^ 2
Local $c = 2 * ASin(Sqrt($a))
Local $d = $R * $c
ConsoleWrite(" 1208 dist - " & $d & @CRLF)
Return $d
EndFunc ;==>_LLDist

If I run script 2 by itself, I get correct console output:

1197 41.38123769313097, -72.50517829321325; 41.38124062679708, -72.50517854467034 1208 dist - 0.000176502556754481

If I run script 1 on the same data, I get:

1197 41.38123769313097, -72.50517829321325; 41.38124062679708, -72.50517854467034 1208 dist - -1.#IND

As you can see, the input data is the same but the output is not.

Additionally, after this happens, other data in some arrays seems to get corrupted, and, for example, an array value that's supposed to range from 1-20,000 or so shows values like -9223372036854762545. A divide by zero causing a memory corruption? Or something else?

Any ideas?

Edited March 19, 2013 by Dana

[Melba23]

Dana,

Interesting.

Have you tried seeing what the other mathematical manipulations used during the function return in the 2 cases as well as the final result? That might pin down the exact line on which the error occurs.

Just one other point - although I do not see why it should affect the result when using the same data - thisline:

Local $a = (Sin($dlat / 2)) ^ 2 + Cos($lat1) * Cos($lat2) * (Sin($dlon / 2)) ^ 2

could do with some more parentheses to reduce its dependence on operator precedence.

M23

[jchd]

That reminds me of ticket #2068

You should use:

Global Const $Pi = 3.1415926535897932385

only once at code init, then use $Pi to keep away from this bug. Perhaps your more complex code (2 UDFs) is invoking ACos(-1) more than 6 times and triggers the issue.

Then I tried a quick'n'dirty Mathematica session to compute the actual result we should get if things were close to perfect:

In[2]:= SetPrecision[{coords={lat1,lon1,lat2,lon2}=(# Pi/180) & /@{41.38123769313097,-72.50517829321325,41.38124062679708,-72.50517854467034},
{dlat=lat2-lat1,dlon=lon2-lon1},
{a=Sin[dlat/2]^2+Cos[lat1] Cos[lat2] Sin[dlon/2]^2,
c=2 ArcSin[Sqrt[a]]},
d=3440.064794816415 c},
20]
Out[2]= {
{41.381237693130970001,-72.505178293213248253,41.381240626797080040,-72.505178544670343399},{0.72223884629551826642,-1.2654540859620937443,0.72223889749765102231,-1.2654540903508479932},
{5.1202132755889806504*10^-8, -4.3887542489073894103*10^-9},
{6.5812556737731055044*10^-16,  5.1307916246026236769*10^-8},
0.00017650255637334405522}

As you can see I requested 20 digits of guaranteed precision for the whole block.

The fact that the points are very close does not seem to impact result accuracy too much, despite here the value of a and c (your $a and $c) are fairly close to 0. Thanks to the "haversine" trick and larger than average precision we don't experience numerical instability in this case, but it might show up silently for closer points or when using standard IEEE floating-point like AutoIt does.

To avoid spurious instability, you can consider points too close as points in a plane instead of a sphere, in which case the error made is unsignificant. I kindly live it to you to compute which conditions make machine-precision rounding error equivalent to sphere-to-plane mapping systematic error.

Also beware of [close to] antipodal points!

Edited March 20, 2013 by jchd

[Spiff59]

I posted a bugtracker a while back describing this:

http://www.autoitscript.com/trac/autoit/ticket/2068

The work-around might be replacing the following:

$lat1 = $lat1 * (ACos(-1) / 180)
$lon1 = $lon1 * (ACos(-1) / 180)
$lat2 = $lat2 * (ACos(-1) / 180)
$lon2 = $lon2 * (ACos(-1) / 180)

With a global at the top of your script like:

$iACos180 = ACos(-1) / 180

And then:

$lat1 = $lat1 * $iACos180
$lon1 = $lon1 * $iACos180
$lat2 = $lat2 * $iACos180
$lon2 = $lon2 * $iACos180

If you've other "ACos(-1)" statements without the division by 180 then maybe just make a global for "ACos(-1)".

But if you execute six "ACos(-1)" commands in row... Blam!

Edit: As Johnny Carson used to say "Weird stuff!"

Edit2: The ticket does claim it was fixed in 3.3.9.0... and testing shows it has been. So leave your code intact and run the beta if you like.

Edited March 20, 2013 by Spiff59

[Dana]

Strange, I hadn't heard of the Acos(-1) bug. I'll check if I'm using the pre or post bugfix version. I had considered starting with $pi = Acos(-1) and then just using $pi thereafter, will try that tomorrow. I am calling the _LLDist function repeatedly, up to maybe a hundred times (I first sort the big database of 20K entries to quickly find points within 1 minute of latitude, usually 60 or so but it varies, then use the function check the actual distance for each of those points, within 1 mile vertially but could be thousands of miles away horizontally, and find the closest).

The odd thing is that the script 2 version uses exactly the same logic, same number of iterations, but doesn't show the problem. Script 1 invokes script 2 only once, after which it's the same.

[jchd]

What you describe seems a good job for the Rtree extension of SQLite.

Don't forget the limited accuracy of FP I mentionned! For really close points, a small perturbation in lattitudes or longitudes can result in a large variation in the (meaningless in that case) result.

Edited March 20, 2013 by jchd

[Dana]

Well, glory be! Globally defining $pi = Acos(-1) at the beginning of the main scrlipt and using $pi in the calculation function fixed everything (except the zillions of ConsoleWrites scattered throughout the file from the debug attempt). Presumably upgrading from 3.3.6.1 might have fixed it, too, but 3.3.9 is still beta (for a year now?).

Still can't see why it behaved differently between the two implementations, unless it were a memory usage thing and I was already on the hairy edge (both of the orignal scripts have a fair amount of array manipulation).

jhcd, this is a brute force approach, but it's adequate for the task. The closest two points that I need to differentiate are some 500' apart, and I don't require an absolute distance, just which point from the database is closed to the GPS derived point. Actually, what takes longest is reading (3 seconds) the [up to] 20,000 line input text file (the file extension is .xls but it's really a .csv file, which I have no control over) and sorting it (_ArraySort, 4 seconds) by latitude for the first rough filter. Using a "top half, bottom half" iterative sort for points within 1 mile north/south and only using the _LLDist function on the 50-100 pre-sort results is quite quick.