[GRASS-dev] [GRASS GIS] #1575: r.horizon bugfix and speedup

[GDAL]

#1575: r.horizon bugfix and speedup
-------------------------+--------------------------------------------------
Reporter: sprice | Owner: grass-dev@…
     Type: enhancement | Status: new
Priority: normal | Milestone: 7.0.0
Component: Default | Version: svn-trunk
Keywords: | Platform: Unspecified
      Cpu: OSX/Intel |
-------------------------+--------------------------------------------------
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

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

--
Ticket URL: <http://trac.osgeo.org/grass/ticket/1575&gt;
GRASS GIS <http://grass.osgeo.org>

[GDAL]

#1575: r.horizon bugfix and speedup
-------------------------+--------------------------------------------------
Reporter: sprice | Owner: grass-dev@…
     Type: enhancement | Status: new
Priority: normal | Milestone: 7.0.0
Component: Default | Version: svn-trunk
Keywords: | Platform: Unspecified
      Cpu: OSX/Intel |
-------------------------+--------------------------------------------------

Comment(by 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;
         }

}}}

--
Ticket URL: <http://trac.osgeo.org/grass/ticket/1575#comment:1&gt;
GRASS GIS <http://grass.osgeo.org>

[GDAL]

#1575: r.horizon bugfix and speedup
-------------------------+--------------------------------------------------
Reporter: sprice | Owner: grass-dev@…
     Type: enhancement | Status: new
Priority: normal | Milestone: 7.0.0
Component: Raster | Version: svn-trunk
Keywords: r.horizon | Platform: All
      Cpu: OSX/Intel |
-------------------------+--------------------------------------------------
Changes (by hamish):

  * keywords: => r.horizon
  * platform: Unspecified => All
  * component: Default => Raster

Old description:

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

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

New description:

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

--

--
Ticket URL: <https://trac.osgeo.org/grass/ticket/1575#comment:2&gt;
GRASS GIS <http://grass.osgeo.org>

[GDAL]

#1575: r.horizon bugfix and speedup
-------------------------+--------------------------------------------------
Reporter: sprice | Owner: grass-dev@…
     Type: enhancement | Status: new
Priority: normal | Milestone: 7.0.0
Component: Raster | Version: svn-trunk
Keywords: r.horizon | Platform: All
      Cpu: All |
-------------------------+--------------------------------------------------
Changes (by neteler):

  * cpu: OSX/Intel => All

Comment:

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

--
Ticket URL: <http://trac.osgeo.org/grass/ticket/1575#comment:3&gt;
GRASS GIS <http://grass.osgeo.org>

[GDAL]

#1575: r.horizon bugfix and speedup
-------------------------+--------------------------------------------------
Reporter: sprice | Owner: grass-dev@…
     Type: enhancement | Status: new
Priority: normal | Milestone: 7.0.0
Component: Raster | Version: svn-trunk
Keywords: r.horizon | Platform: All
      Cpu: All |
-------------------------+--------------------------------------------------

Comment(by 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!

--
Ticket URL: <https://trac.osgeo.org/grass/ticket/1575#comment:4&gt;
GRASS GIS <http://grass.osgeo.org>

[GDAL]

#1575: r.horizon bugfix and speedup
-------------------------+--------------------------------------------------
Reporter: sprice | Owner: grass-dev@…
     Type: enhancement | Status: new
Priority: normal | Milestone: 7.0.0
Component: Raster | Version: svn-trunk
Keywords: r.horizon | Platform: All
      Cpu: All |
-------------------------+--------------------------------------------------

Comment(by 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.

--
Ticket URL: <https://trac.osgeo.org/grass/ticket/1575#comment:5&gt;
GRASS GIS <http://grass.osgeo.org>

[GDAL]

#1575: r.horizon bugfix and speedup
-------------------------+--------------------------------------------------
Reporter: sprice | Owner: grass-dev@…
     Type: enhancement | Status: new
Priority: normal | Milestone: 7.0.0
Component: Raster | Version: svn-trunk
Keywords: r.horizon | Platform: All
      Cpu: All |
-------------------------+--------------------------------------------------

Comment(by 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

--
Ticket URL: <http://trac.osgeo.org/grass/ticket/1575#comment:6&gt;
GRASS GIS <http://grass.osgeo.org>