r.horizon bugfix and speedup

[sprice]

I mentioned this on the mailing list, but didn't get a response. So here's a more official report:

I made some tweaks to r.horizon/main.c to almost half its runtime and remove a bug. By my tests, it's output is identical to before, but a 219 minute run has been reduced to 110 minutes. Mostly due to removal of unneeded floor() calls. I could use some people to test and make sure I'm not introducing any bugs, though. Thanks, Seth

[sprice]

Well let's try that patch again...

I made some tweaks to r.horizon/main.c to almost half its runtime and
remove a bug. By my tests, it's output is identical to before, but a 219
minute run has been reduced to 110 minutes. Mostly due to removal of
unneeded floor() calls. I could use some people to test and make sure I'm
not introducing any bugs, though.
Thanks,
Seth

--- main.orig	2011-08-14 06:32:49.000000000 -0600
+++ main.c	2012-02-12 19:20:15.000000000 -0700
@@ -144,12 +144,15 @@
/* use G_distance() instead ??!?! */
double distance(double x1, double x2, double y1, double y2)
{
+	double diffX = x1 - x2;
+	double diffY = y1 - y2;
+
    if (ll_correction) {
-	return DEGREEINMETERS * sqrt(coslatsq * (x1 - x2) * (x1 - x2)
-				     + (y1 - y2) * (y1 - y2));
+	return DEGREEINMETERS * sqrt(coslatsq * diffX * diffX
+				     + diffY * diffY);
    }
    else {
-	return sqrt((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2));
+	return sqrt(diffX * diffX + diffY * diffY);
    }
}

@@ -841,36 +844,28 @@

int new_point()
{
-    int iold, jold;
-    int succes = 1, succes2 = 1;
-    double sx, sy;
-    double dx, dy;
-
-    iold = ip;
-    jold = jp;
+    int iold = ip;
+    int jold = jp;

-    while (succes) {
+    while (1) {
	yy0 += stepsinangle;
	xx0 += stepcosangle;


	/* offset 0.5 cell size to get the right cell i, j */
-	sx = xx0 * invstepx + offsetx;
-	sy = yy0 * invstepy + offsety;
-	ip = (int)sx;
-	jp = (int)sy;
+	ip = (int) (xx0 * invstepx + offsetx);
+	jp = (int) (yy0 * invstepy + offsety);

	/* test outside of raster */
	if ((ip < 0) || (ip >= n) || (jp < 0) || (jp >= m))
	    return (3);

	if ((ip != iold) || (jp != jold)) {
-	    dx = (double)ip *stepx;
-	    dy = (double)jp *stepy;
+	    double dx = (double)ip *stepx;
+	    double dy = (double)jp *stepy;

	    length = distance(xg0, dx, yg0, dy);  /* dist from orig. grid point
to the current grid point */
-	    succes2 = test_low_res();
-	    if (succes2 == 1) {
+	    if (test_low_res() == 1) {
		zp = z[jp][ip];
		return (1);
	    }
@@ -882,54 +877,44 @@

int test_low_res()
{
-    int iold100, jold100;
-    double sx, sy;
-    int delx, dely, mindel;
-    double zp100, z2, curvature_diff;
-
-    iold100 = ip100;
-    jold100 = jp100;
-    ip100 = floor(ip / 100.);
-    jp100 = floor(jp / 100.);
+    int iold100 = ip100;
+    int jold100 = jp100;
+    ip100 = ip * .01;
+    jp100 = jp * .01;
    /*test the new position with low resolution */
    if ((ip100 != iold100) || (jp100 != jold100)) {
	/* G_debug(3,"%d %d %d %d\n",ip,jp, iold100,jold100); */
	/*  replace with approximate version
	   curvature_diff = EARTHRADIUS*(1.-cos(length/EARTHRADIUS));
-	 */
+
+	Code folded into the 'if' statement:
	curvature_diff = 0.5 * length * length * invEarth;
	z2 = z_orig + curvature_diff + length * tanh0;
	zp100 = z100[jp100][ip100];
+	if (zp100 <= z2) {...
+	 */
	/*G_debug(3,"%d %d %lf %lf \n",ip,jp,z2,zp100); */

-	if (zp100 <= z2)
+	if (z100[jp100][ip100] <= z_orig + 0.5 * length * length * invEarth +
length * tanh0)
	    /*skip to the next lowres cell */
	{
-	    delx = 32000;
-	    dely = 32000;
+	    int mindel;
+	    int delx = 32000;
+	    int dely = 32000;
+		double sx = (xx0 * invstepx + offsetx) * .01;
+		double sy = (yy0 * invstepy + offsety) * .01;
+
	    if (cosangle > 0.) {
-		sx = xx0 * invstepx + offsetx;
-		delx =
-		    floor(fabs
-			  ((ceil(sx / 100.) - (sx / 100.)) * distcosangle));
+		delx = fabs((ceil(sx) - sx) * distcosangle);
	    }
-	    if (cosangle < 0.) {
-		sx = xx0 * invstepx + offsetx;
-		delx =
-		    floor(fabs
-			  ((floor(sx / 100.) - (sx / 100.)) * distcosangle));
+	    else if (cosangle < 0.) {
+		delx = fabs((floor(sx) - sx) * distcosangle);
	    }
	    if (sinangle > 0.) {
-		sy = yy0 * invstepy + offsety;
-		dely =
-		    floor(fabs
-			  ((ceil(sy / 100.) - (sy / 100.)) * distsinangle));
+		dely = fabs((ceil(sy) - sy) * distsinangle);
	    }
	    else if (sinangle < 0.) {
-		sy = yy0 * invstepy + offsety;
-		dely =
-		    floor(fabs
-			  ((floor(jp / 100.) - (sy / 100.)) * distsinangle));
+		dely = fabs((floor(sy) - sy) * distsinangle);
	    }

	    mindel = min(delx, dely);
@@ -953,17 +938,14 @@

double searching()
{
-    double z2;
-    double curvature_diff;
-    int succes = 1;
-
    if (zp == UNDEFZ)
	return 0;

    while (1) {
-	succes = new_point();
+    double z2;
+    double curvature_diff;

-	if (succes != 1) {
+	if (new_point() != 1) {
	    break;
	}

[neteler]

I have manually merged the inline diff line by line (gngn) and made it an attachment.

[hamish]

Hi,

C Q: is declaring variables mid-function (say at the start of an if(){} statement) fully portable to all the platforms/compilers that we wish to support? how about declaring and initializing variables at the same time? +/-'s?

does the /100 -> *0.01 remove the expensive divide op, or is it just about removing the cumulative number of ops?

is declaring a new variable relatively inexpensive vs declaring it then using it?

which part of the patch specifically is the bug?

(as concerned with readability/maintainability as I am with speed; w/ code comments where appropriate we can have the best of both worlds; the "Code folded into the 'if' statement:" explanation is great tx)

thanks, educational-mode non-comp.sci. major Hamish

ps- in general I think we can gain a huge amount from running various bits of grass through a profiler. the "Bugs" wiki page has some suggestions, the yes,still-free-for-non-commercial/research-use Intel C compiler has a very good profiler+tutorial too. thanks Seth!

[sprice]

As far as I know, declaring variables in the if(){} is portable. I'm not sure what the speed advantage is, though. I like to explicitly give the compiler the opportunity to optimize, though.

The *0.01 is to remove the divide op.

I'm not sure where the bug is. It's been a while.

[neteler]

patch updated to relbranch7 r64094

[neteler]

patch updated to relbranch7 r64094

[wenzeslaus]

Probably a thorough test of current and proposed r.horizon is needed before somebody will change the existing implementation. Please see: ​http://grass.osgeo.org/grass71/manuals/libpython/gunittest_testing.html

[martinl]

Milestone renamed

[martinl]

Seems to be solved, closing.