Opened 3 years ago

Last modified 2 months ago

#3587 new defect

lwgeom_mindistance2d slows down topology loading

Reported by: strk Owned by: strk
Priority: medium Milestone: PostGIS Fund Me
Component: topology Version: 2.2.x
Keywords: performance Cc: dbaston

Description

Profiling an slow addition of a single point I found the CPU busy in computing the minimum distance between the newly added edge with 43520 vertices and each of 2862 other edges intersecting it. Such operation takes a very long time, even multiple seconds for each single distance computation.

Overall, the addition of a that _single_ line takes over 13 minutes.

The same long timings can be seen using ST_Distance itself. Examples:

LOG:  duration: 1043.816 ms  statement: select a.id aid, b.id bid, st_npoints(a.wkb_geometry), st_npoints(b.wkb_geometry), st_distance(a.wkb_geometry, b.wkb_geometry) from contour_onet a, contour_onet b where a.id = 14408 and b.id = 9928;
  aid  | bid  | st_npoints | st_npoints |   st_distance    
-------+------+------------+------------+------------------
 14408 | 9928 |      43520 |        499 | 10568.8576870335
(1 row)
LOG:  duration: 4937.063 ms  statement: select a.id aid, b.id bid, st_npoints(a.wkb_geometry), st_npoints(b.wkb_geometry), st_distance(a.wkb_geometry, b.wkb_geometry) from contour_onet a, contour_onet b where a.id = 14408 and st_npoints(b.wkb_geometry) = 2247;
  aid  | bid  | st_npoints | st_npoints |   st_distance
-------+------+------------+------------+------------------
 14408 | 9031 |      43520 |       2247 | 53917.1247457271
(1 row)

This ticket is to find speed improvement opportunities.

Attachments (4)

t.sql.tar.bz2 (1.5 MB) - added by nicklas 3 years ago.
bypassBboxCheckAndFastDistanceCalc.diff (1.2 KB) - added by nicklas 3 years ago.
A patch that comments out the code that performs the bbox test and sends to faster distance calc
slodist.sql.gz (493.1 KB) - added by strk 3 years ago.
A dump of the "slowdist" subject table
slowdist.png (20.4 KB) - added by strk 3 years ago.

Change History (42)

comment:1 Changed 3 years ago by strk

Another example:

LOG:  duration: 6416.539 ms  statement: select st_summary(a), st_summary(b), st_distance(a,b) from slowdist;
            st_summary            |           st_summary            |   st_distance   
----------------------------------+---------------------------------+-----------------
 LineString[BS] with 43520 points | LineString[BS] with 3153 points | 254918.45088224
(1 row)

comment:2 Changed 3 years ago by strk

In 14977:

Make adding a line to topology interruptible

... as it otherwise gets frustrating when slow (see #3587)

Changed 3 years ago by nicklas

Attachment: t.sql.tar.bz2 added

comment:3 Changed 3 years ago by nicklas

Yes, this is good.

I think Paul have the solution, we just have to encourage him a little.

Some background: My first contribution to PostGIS was about this, added in 1.5.

But the result of that only helps when bounding boxes is not overlapping. What I did is described in: https://trac.osgeo.org/postgis/wiki/NewDistCalcGeom2Geom

This gives a good boost, when bbox isn't overlapping.

I added 3 geometries as example in attachement.

I copied and moved one of them so the bounding boxes overlap.

So those 2 queries is using the same geometries, but the first query with overlapping bounding boxes:

select st_npoints(a.geom), st_npoints(b.geom), st_distance(a.geom, b.geom) from 
(select geom from t where gid = 1) a,
(select geom from t where gid = 2) b

st_npoints | st_npoints |   st_distance   
------------+------------+-----------------
      55773 |      31248 | 1791.9837591391
(1 row)


Time: 73239.353 ms
select st_npoints(a.geom), st_npoints(b.geom), st_distance(a.geom, b.geom) from 
(select geom from t where gid = 1) a,
(select geom from t where gid = 3) b
 st_npoints | st_npoints |   st_distance    
------------+------------+------------------
      55773 |      31248 | 194387.282408659
(1 row)

Time: 709.125 ms

So this is good, but useless in many situations. But Paul have some tree-structure that will do the thing.

Do it Paul, just do it :-)

comment:4 Changed 3 years ago by strk

Is there any way to access the older algorithm, for test ? I think the tree structure you are talking about is the one in postgis/ ? Or is that one for the new algorithm you are talking about ? Whatever it is, having access to that tree in liblwgeom would greatly help...

comment:5 Changed 3 years ago by nicklas

That tree thing, I think Paul have implemented in geography, but I am no sure if it is the same thing he talked about to implement in geometry 5 years ago.

About the old algorithm, I don't get what you mean. In 1.4 and before all distance calculations was done eith a "brute force" method. Just comparing all combinations of vertex-vertex and vertex-edge. In 1.5 I implemented a more efficient method for cases with linestring and polygons involved, and bboxes that don't intersect. For cases with overlapping bboxes it falls back to the same old "brute force" method, so it is still there. Not exectly the same as in 1.4, since I did some smaller changes. Since 1.5 I think the distance code for geometry is more or less is untouched except that someone (Paul I think) have added support for circular geometries.

comment:6 Changed 3 years ago by strk

In my case I'm only calling lwgeom_mindistance2d against linestring pairs with intersecting bounding boxes. So you're saying I do am seeing the brute force in action ? No other overhead other than a re-check for bounding box overlap ?

Paul: is the rtree cache used for ST_Distance(geometry, geometry) ? Any plan to move the rtree interface down to liblwgeom for use from the topology code ? (it is common to test the same long line against hundreds of other lines, so a local temporary cache could be very useful)

comment:7 Changed 3 years ago by pramsey

There is no rtree cache used for distance in geometry currently. The start of the code is in liblwgeom/lwtree.c but it has never been hooked up and is probably rotten by now. I got pulled into circular geometry distance by the requirements of the tree, and never got back to it. The tree code in geography turned up enough cases where brute_distance() != tree_distance() that I have viewed hooking up the tree code as a huge task requiring piles of testing and and QA/QC. It would certainly be very useful for repeated distance testing, where the repeated side is complex.

The only place where we have tree code that has not been pulled down into liblwgeom is the tree code for point-in-poly intersects true/false testing.

comment:8 Changed 3 years ago by nicklas

Yes, you are meeting the brute force wall at once the bboxes is overlapping.

About the overhead, things can probably be optimized. But my guess is that will not give more than maybe 10 % increase of speed. If I recall right I had the brute force part run approx 10 % faster when I rewrote it for handling distance and shortest line the same way. The purpose was to get exactly the same distance as the length of the shortest line. As a bonus there was a small increase in speed.

comment:9 Changed 3 years ago by strk

I'm still surprised by the excess time. Note that using ST_Intersects(ST_Buffer(a,10), b) takes about 200ms. Given the topology use for mindistance2d is to find out if the geometry is within a given amount of units from each other, I guess the check could be replaced with those two calls to GEOS.

But the extreme difference in time seems to suggest there's something fundamentally wrong with the distance algorithm. Is there any way to access the older algorithm, fully brute force, from SQL ?

comment:10 Changed 3 years ago by strk

See the numbers:

LOG:  duration: 6360.923 ms  statement: select st_dwithin(a,b,10) from slowdist;
 st_dwithin 
------------
 f
LOG:  duration: 210.333 ms  statement: select st_intersects(st_buffer(a,10),b) from slowdist;
 st_intersects 
---------------
 f
LOG:  duration: 59.426 ms  statement: select st_intersects(a,st_buffer(b,10)) from slowdist;
 st_intersects 
---------------
 f

comment:11 Changed 3 years ago by nicklas

Yes, I agree it looks very strange. If you want to test the version before my contribution you should fire up a 1.4.

It should be fairly easy to force the brute force method in todays code too. I can take a look where to search. But that will not rewind to what it was before I added things.

I will try to take a look too if I can see what is happening.

Can the prepared geometry help the intersects case, or is it only when quering many geometries that help?

comment:12 Changed 3 years ago by strk

The numbers you see are w/out prepared geometries, I think, because a single ST_Intersects call is issued within the statement.

comment:13 Changed 3 years ago by strk

PostGIS 1.4 is actually worst:

pg14=# select st_intersects(a,st_buffer(b,10)) from slowdist;
 st_intersects 
---------------
 f
(1 row)

Time: 69.658 ms
pg14=# select st_intersects(st_buffer(a,10),b) from slowdist;
 st_intersects 
---------------
 f
(1 row)

Time: 423.872 ms
pg14=# select st_dwithin(a,b,10) from slowdist;
 st_dwithin 
------------
 f
(1 row)

Time: 13152.853 ms

Changed 3 years ago by nicklas

A patch that comments out the code that performs the bbox test and sends to faster distance calc

comment:14 Changed 3 years ago by nicklas

I have attached a patch that comments out the bbox test and sending to faster distance calc.

I see no difference in timing with that.

I also think this is very strange.

I can see 3 possible reasons

1) Something very smart is happening in GEOS in st_intersects code

2) We have some very big overhead in our helping functions that for instance do very many memory allocations or reallocations, that I don't realize when I look at the code.

3) There is some ugly brain fart in the distance code itself. But since 1.4 shows the same pattern, I guess that is something we have had for a very long time and that I failed to see when I did my changes.

An interesting thing is that ST_3DIntersects behaves in about the same way as the distance functions. That is not strange since it is a wrapper around similar code as the 2D distance functions.

comment:15 Changed 3 years ago by pramsey

My bet is on (1), as I recall a conversation w/ Martin in which is said his distance implementation used a sweepline implementation. He may be tossing all the edges into a bin and discarding many of them far more efficiently than we are. Most of your enhancements were around making multipart objects more efficient, but I think strk's examples are basically one huge part, right? It's possible my tree implementation would be as good as Martin's implementation, but it might not necessarily, as a sweepline could get a more consistent result than a depth-first search, which I'm pretty sure is what I'm getting out of my recursing on the trees.

Changed 3 years ago by strk

Attachment: slodist.sql.gz added

A dump of the "slowdist" subject table

comment:16 Changed 3 years ago by strk

I've attached the "slodist" table used for my times above, so you can reproduce. I also bet on GEOS being smart, but note this is a relate operation, not a distance operation. You get the same fast times with ST_Relate(ST_Buffer(..))

comment:17 Changed 3 years ago by strk

I confirm objects are single-part:

comment:18 Changed 3 years ago by nicklas

small correction @pramsey, the faster way of doing it is all about single geometries. But I couldn't fins a way of using it when geometries is too close.

Interesting also that when I tried it on your geoemtries @strk, it really showed that it doesn't do a goog work when such a long geometry as a is close to the smaller one.

I tried by translating the geometry away from the bounding box to make the faster algorithm kick in. When I move it just outsinde the bounding box, the timing was even worse than the brute force version. But when I move it further away it start to kick ass.

without moving anything on my box:

slowdist=# select st_distance(a, b) from slowdist;
   st_distance   
-----------------
 254918.45088224
(1 row)

Time: 7134.470 ms

Just outside the bbox (worst case inside the "circle" of the long line):

slowdist=# select st_distance(a, st_translate(b,-290000,0)) from slowdist;
   st_distance    
------------------
 432666.185703918
(1 row)


Time: 9459.979 ms

But if I move the geometries more apart so geometry a is more in one direction from b (happens easier with geoemtries that is not so long as a) then we see the difference:

slowdist=# select st_distance(a, st_translate(b,-29000000,0)) from slowdist;
   st_distance    
------------------
 28728819.1262113
(1 row)

Time: 23.064 ms

I haven't seen before cases when the "faster" way is slower than brute force, but I guess this is a little extreme case.

If we insted move geometry b just above the bounding box of a we get nicer results at once:

slowdist=# select st_distance(a, st_translate(b,0,1950000)) from slowdist;
   st_distance    
------------------
 79242.5337417284
(1 row)

Time: 250.828 ms

This doesn't help the case.

But if @pramsey means that geos have also distance algorithms with the performance of ST_Intersects, I guess we should switch to them at once.

If geos (jts) is that good also at distance calculations I think that is something Martin has implemented in the last years. I remember Martin reviewing my work and mentioned that he had thought in the same direction but not tested it. He also confirmed that it will not work on too close geometeries.

comment:19 Changed 3 years ago by nicklas

In strks example we get 43520 * 3153 combinations to test. For all those combinations we have to test both point to point and point to segment.

So I guess it is possible that our computers isn't faster than this. If using 7 seconds for the distance caclulation it means about 19000 combinations per millisecond. 19000 point to point and 19000 point to segment (well I guess it is two less point to segment combination).

But there might be things to do about it also with brute force. I mean if we can decrease with 1 microsecond per combination it will give a huge effect.

When doing it the brute force way, I think it is 2 places to look for optimization.

1) Both for point/point and point/segment tests we end up with a sqrt calculation. When I worked on those things long time ago I tested to not do that sqrt calc and only do it on the resulting closest points. But I ran into many strange memory problems. I didn't investigate to the ens, but I guess it results in handling very big numbers which double precision can't handle. I mean 10000 km results in handling the number 10 raised by 14. But if that is solvable (by me just doing something wrong back then) I guess that could make a noticable difference to avoid 38000 sqrt calculations per millisecond.

2) at line 2075 in measures.c is a calculation with a lot of addition, substraction, multiplication and division. That line is executed for each point/segment test. There is more below in the same function, but that is only executed when the projected point on the segment is actually on the segment. More common is that the projected point is outside the segment, and the function goes directly to compare the endpoint of the segment with the point.

comment:20 Changed 3 years ago by strk

There's another tree implementation in lwgeom_geos_cluster.c, simply wrapping the GEOS STRTree class. In case we want to implement something like that internally. Or GEOS C-API does have a GEOSDistance function if we want to try that.

About avoiding the sqrt, the geometries could be shifted to the origin to lower the numbers

comment:21 Changed 3 years ago by strk

FYI: the GEOSDistance implementation is slower than the PostGIS one, for this case: 30 seconds vs. 15 seconds (running on battery).

For sure it's not the same code path issued by ST_Intersects...

Changed 3 years ago by strk

Attachment: slowdist.png added

comment:22 Changed 3 years ago by strk

The offending lines: A is the long one, B the small one

comment:23 Changed 3 years ago by strk

I confirm the MonotoneChain? / Sweepline based algorithm is used for constructing the lines topologies (for ST_Intersects / ST_Relate).

If we take the STRTree route, GEOS-3.6.0 (unreleased) will have a GEOSSTRtree_nearest function, so we could build a tree containing the envelopes of each segment (or monotonechain) of the longest geometry and then query it for closest segment (or chain) on each of the shorter geometry segments.

comment:24 Changed 3 years ago by pramsey

No, you don't want to build an strtree of the edges, as that will take quite a bit of computational time itself. The trick to trees is to use the natural autocorrelation of the edges to group the leaves into internal nodes. That way no sorting or individual insertion is required, you just roll up the ordered set of edges into a tree in O(N) time, then you recurse back down the tree to find the answer in O(log(N)) time, so in theory the whole thing is never more than O(N) + O(log(N)) = O(N). Now, why Intersects() is still quite fast remains a mystery.

comment:25 Changed 3 years ago by strk

I've just tried creating the STRtree and it does not seem to affect the timings much. No querying, just creation and insert. The time is still ~15 seconds.

comment:26 Changed 3 years ago by strk

querying the STRtree for nearest segment gets to the distance result in ~20ms. Looks promosing. If only the result was similar to the original one ...

LOG:  duration: 19.001 ms  statement: select st_summary(a), st_summary(b), st_distance(a,b) from slowdist;
            st_summary            |           st_summary            |   st_distance
----------------------------------+---------------------------------+-----------------
 LineString[BS] with 43520 points | LineString[BS] with 3153 points | 827058.14279572
(1 row)

comment:27 Changed 3 years ago by strk

Ok, a correct implementation takes ~2.3 seconds. It's not 20ms but still much faster than ~13 seconds...

LOG:  duration: 2320.946 ms  statement: select st_distance(b,a) from slowdist;
   st_distance   
-----------------
 254918.45088224
(1 row)
LOG:  duration: 2330.856 ms  statement: select st_dwithin(b,a,10) from slowdist;
 st_dwithin 
------------
 f
(1 row)

It just needs to be cleaned up to go.

Pull request follows:

The following changes since commit f4d8771b218fcf2fc66db3f1ffeecda3b4e8a140:

  Remove trailing spaces (2016-07-04 23:50:26 +0200)

are available in the git repository at:

  https://git.osgeo.org/gogs/postgis/postgis.git svn-trunk-strtree-mindist

for you to fetch changes up to c1e163fa3883b75094623cdd66335f9d4d973727:

  Use STRtree based mindistance computation (2016-07-05 01:12:20 +0200)

----------------------------------------------------------------
Sandro Santilli (1):
      Use STRtree based mindistance computation

 liblwgeom/measures.c | 62 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++-
 1 file changed, 61 insertions(+), 1 deletion(-)

Tests are welcome!

Regina, can we find a way to point Jenkins at Gogs pull requests, maybe ?

comment:28 Changed 3 years ago by strk

Cc: dbaston added

Daniel, given I'm getting a segment which doesn't seem to be the closest one from GEOSSTRtree_nearest, maybe you want to take a look at the PR too ? There could be a bug in GEOS.

What I'm doing:

  1. Creating an STRTree with node capacity 10
  2. Populating the tree with 43519 segments
  3. Finding the closest segment to each of 3152 segments of the other geometry
  4. Computing the closest points within the segment pair

The problem is that the brute force approach finds a pair of segments which are closer than the ones found by GEOSSTRtree_nearest

comment:29 Changed 3 years ago by strk

Daniel: forget it, it was an error in my code, now works correctly, ad it is still 15 times faster than before:

LOG:  duration: 954.569 ms  statement: select ST_AsText(st_shortestline(a,b)) from slowdist;
                              st_astext                              
---------------------------------------------------------------------
 LINESTRING(5161330.6 -4881460.77503163,5046321.13099501 -4653960.8)
(1 row)

Not as good as ST_Intersection, but still better than what we have now. I guess some experimenting could be done with different nodes capacity (10 sounds a bit low)

comment:30 Changed 3 years ago by strk

An STRtree node capacity of 256 resulted in raising the timings to ~5 seconds. A node capacity of 32 made timings =~ 2 seconds. A node capacity of 4 made timings =~ 686 ms. A node capacity of 2 made timings =~ 532 ms. A node capacity of 1 resulted in a segfault.

I'm not sure what to make of these numbers.

comment:31 Changed 3 years ago by dbaston

If you're willing to require GEOS 3.6, we should probably port IndexedFacetDistance? over from JTS. From a quick glance, it looks like all of the non-trivial code was already captured with the MinimumClearance? port. This would avoid needing to build a GEOS LineString? for every line segment.

comment:32 Changed 3 years ago by strk

Good point about the LineString?-per-segment need removal. I've to say IndexedFacetDistance? also already has a good API for use from the topology code, so good idea !

comment:33 Changed 3 years ago by strk

Another great idea (for topology): DROP the distance computation. Given that distance is only used to reduce the number of ST_Snap operations, and that ST_Snap would be a no-op if no vertex is closer than the given tolerated distance... we might as well just avoid to perform the pre-flight check.

The testsuite is fine with that, which is a good sign!

comment:34 Changed 3 years ago by dbaston

FWIW, I ported IndexedFacetDistance? to GEOS and exposed it in the CAPI as GEOSDistanceIndexed. Using the two geometries from the table strk posted, the calculation timings are as follows (includes overhead of reading text file and parsing WKT)

GEOSDistance: 41.7 s

GEOSDistanceIndexed: 0.084s

This might not be the right fit for the situation described in this ticket, but it seems like it could be useful to PostGIS overall.

comment:35 Changed 3 years ago by robe

Milestone: PostGIS 2.2.3PostGIS 2.4.0

comment:36 Changed 2 years ago by robe

Milestone: PostGIS 2.4.0PostGIS 2.5.0

comment:37 Changed 15 months ago by robe

Milestone: PostGIS 2.5.0PostGIS 3.0.0

comment:38 Changed 2 months ago by komzpa

Milestone: PostGIS 3.0.0PostGIS Fund Me
Note: See TracTickets for help on using tickets.