Transcript Flow Computation on Massive Grids

External Memory Graph Algorithms
and
Applications to GIS
Laura Toma
Bowdoin College
Outline
 A GIS application: TerraFlow
 I/O-efficient graph algorithms
• Problems and results
• Graph separation
• Topological sort on planar DAGs
 Future directions/open questions
Why GIS?
 How it all started..
• Duke environmental researchers
• Computing flow accumulation for Appalachian Mountains took 14
days (with 512MB memory)
• 800km x 800km at 100m resolution  ~64 million points
 GIS (Geographic Information Systems)
• System that handles spatial data
• Visualization, processing, queries, analysis…
• Rich source of problems for Computer Science
• Graphics, graph theory, computational geometry, scientific
computing…
GIS and the Environment
Precipitation in Tropical
South America
Indispensable tool
•
Monitoring: keep an eye on the
state of earth systems using
satellites and monitoring stations
(water, pollution, ecosystems,
urban development,…)
•
Modeling and simulation: predict
consequences of human actions
and natural processes
• Analysis and risk assessment: find
the problem areas and analyse the
possible causes (soil erosion,
groundwater pollution,..)
•
Planning and decision support:
provide information and tools for
better management of resources
Lots of rain
Dry
Nitrogen in Chesapeake Bay
High nitrogen
concentrations
Computations on Terrains
Reality:
Elevation of terrain is a
continuous function of two
variables h(x,y)
Estimate, predict, simulate
 Flooding, pollution
 Erosion, deposition
 Vegetation structure
 ….
GIS:
DEM (Digital Elevation Model) is a
set of sample points and their
heights
{ x, y, hxy}
Model and
compute indices
TerraFlow
DEM Representations
TIN
3
7
7
3
7
7
2
5
1
2
5
1
4
8
9
4
8
9
3
7
7
3
7
7
2
5
1
2
5
1
4
8
9
4
8
9
Grid
Contour lines
Sample points
TerraFlow
Modeling Flow on Terrains
 What happens when it rains?
• Predict areas susceptible to floods.
• Predict location of streams.
 Flow is modeled by computing two basic attributes from the DEM
of the terrain:
• Flow Direction (FD)
• The direction water flows at a point
• Flow Accumulation (FA)
• Total amount of water that flows through a point
if water is distributed according to the flow directions
TerraFlow
Flow Accumulation
of Panama
TerraFlow
Panama Flow Accumulation: zoom
Uses
Flow direction and flow accumulation are used for:

Computing other hydrological attributes
• river network
• moisture indices
• watersheds and watershed divides

Analysis and prediction of sediment
and pollutant movement in landscapes.

Decision support in land management, flood and pollution prevention and
disaster management
Massive Terrain Data
 Remote sensing technology
• Massive amounts of terrain
data
• Higher resolutions (1km,
100m, 30m, 10m, 1m,…)



NASA-SRTM
• Mission launched in 2001
• Acquired data for 80% of earth
at 30m resolution
• 5TB
USGS
• Most of US at 10m resolution
LIDAR
• 1m res
TerraFlow
Example: LIDAR Terrain Data
 Massive (irregular) point sets (1-10m resolution)
 Relatively cheap and easy to collect
Example: Jockey’s ridge (NC coast)
It’s Growing!
 Appalachian Mountains
Area if approx. 800 km x 800 km
Sampled at:
• 100m resolution:  64 million points (128MB)
• 30m resolution:  640
(1.2GB)
• 10m resolution:  6400 = 6.4 billion (12GB)
• 1m resolution:  600.4 billion
(1.2TB)
I/O Model [AV’88]
Block I/O
M
 I/O-operation:
• movement of one block of data from/to disk
 Complexity measure: number of I/Os
 Fundamental bounds:
Scanning: scan(N) = O ( NB )
N = problem size
B = disk block size
M = memory size
Sorting:
sort(N) = O( NB log M B
 In practice B and M are big
N
B
 NB log M B
N
B
 N
I/Os
N
) I/Os
B
Flow Direction (FD) on Grids
 Water flows downhill
• follows the gradient
 On grids: Approximated using 3x3 neighborhood
• SFD (Single-Flow Direction):
• FD points to the steepest downslope neighbor
• MFD (Multiple-Flow direction) :
• FD points to all downslope neighbors
Computing FD




Goal: compute FD for every cell in the grid (FD grid)
Algorithm:
• For each cell compute SFD/MFD by inspecting 8 neighbor cells
Analysis: O(N) time for a grid of N cells
Is this all?
• NO! flat areas: Plateas and sinks
TerraFlow
FD on Flat Areas
 Plateaus
• A cell flows towards the nearest spill point on the boundary of the plateau
• Compute FD on plateaus using CC and BFS
 Sinks
• Route the water uphill out of the sink by modeling
flooding: uniformly pouring water on terrain
until steady-state is reached
• Flooding removes (fills) sinks
 Assign uphill flow directions on the original terrain by assigning
downhill flow directions on the flooded terrain
TerraFlow
Flooding
 Watershed: part of the terrain that flows into a sink
 Sinks  partition of terrain into watersheds  watershed graph GT
• Vertices are watersheds; add vertex for the “outside” watershed
• Edge (u,v) if watersheds u,v are adjacent
• Edge (u,v) labeled with lowest height on boundary between u and v
 Flooding: Compute for each watershed u to the height hu of the lowestheight path in GT from u to the “outside” watershed.
• the height of a path is the height of the highest edge on path
TerraFlow
Flooding
 Plane-sweep algorithm with a Union-Find structure
• Initially only the outside watershed is done
• Sweep watershed graph bottom-up with a horizontal plane
• When hit edge (u,v)
• If both watersheds u and v are done, ignore
• If none is done, union them
• If precisely one is not done, raise it at h(u,v) and mark it done
 Theorem:
Flooding and the FD grid can be computed in O(sort(N)) I/Os on a grid
DEM of size N.
Flow Accumulation (FA) on Grids
FA models water flow through each cell with “uniform rain”
• Initially one unit of water in each cell
• Water distributed from each cell to neighbors pointed to by its FD
• Flow conservation: If several FD, distribute proportionally to
height difference
• Flow accumulation of cell is total flow through it
Goal: compute FA for every cell in the grid (FA grid)
Computing FA
 FD graph
• node for each cell
• edge from a to b if FD of a points to b
 FD graph must be acyclic
 FD graph
• SFD: a set of trees
• MFD: a DAG
TerraFlow
Standard FA Algorithm
 Algorithm:
• Input: flow direction grid FD
• Output: flow accumulation grid FA (initialized to 1)
• Process (sweep) cells in topological order. For each cell:
• Read flow from FA grid and direction from FD grid
• Update flow in FA grid for downslope neighbors
 Analysis
• One sweep enough: O(sort) + O(N) time for a grid of N cells,
• ..but O(N) I/Os: Cells in topological order distributed over the terrain
Scalability Problem
 We can compute FD and FA using simple O(N)-time algorithms
 ..but.. for large sets..??
Dataset Size (log)
I/O-Efficient Flow Accumulation
[ATV00]
 Eliminating scattered accesses to FD grid
• Store FD grid in topological order
 Eliminating scattered accesses to FA grid
• Obs: flow to neighbor cell only needed at the time when the
neighbor is processed:
• Time when cell is processed Topological rank
priority
• Push flow by inserting flow increment in priority queue with
priority equal to neighbor’s time
• Flow of cell obtained using DeleteMin operations
 Use I/O-efficient priority queue [A95, BK97]
N
N
• O(N) operations in O( B log M B B ) I/Os
TerraFlow
http://www.cs.duke.edu/geo*/terraflow
TerraFlow
 TerraFlow: implementation of I/O-efficient FD and FA algorithms
• Significantly faster on very large grids than existing GIS software
• Scalable: 1 billion elements!! (>2GB data)
 Implementation
• C++, uses TPIE (Transparent Parallel I/O Environment)
• Library of I/O-efficient modules developed at Duke
 Experimental platform
• TerraFlow, ArcInfo: 500MHz Alpha, FreeBSD 4.0, 1GB RAM
• GRASS/TARDEM: 500MHz Intel PIII, FreeBSD/Windows, 1GB RAM
http://www.cs.duke.edu/geo*/terraflow
TerraFlow
Experimental Results
 GRASS cannot handle Hawaii dataset (killed after 17 days)
 TARDEM cannot handle Cumberlands dataset (killed after 20 days)
 Significant speedup over ArcInfo (ESRI) for large datasets
90
• East-Coast
TerraFlow: 8.7 Hours
ArcInfo: 78 Hours
TerraFlow 512
80
TerraFlow 128
ArcInfo 512
70
ArcInfo 128
50
ArcInfo
40
TerraFlow
30
20
10
to
n
t
W
as
hi
ng
es
id
w
t-C
Ea
s
M
oa
s
t
E
N
er
rla
be
um
C
Lo
w
nd
s
ai
i
aw
ev
N
er
ra
Si
H
ad
a
o
Ri
c
er
to
Pu
we
a
h
0
Ka
• Washington
TerraFlow: 63 Hours
ArcInfo: %
Running Time (Hours)
60
http://grass.itc.it
TerraFlow
TerraFlow in GRASS
PIII dual 1GHz processor, 1GB RAM
Grid size
r.terraflow
r.watershed
1.85 min
9.2 min
Dataset
Grid dimensions
(million
elements)
Kaweah
1163 x 1424
1.6
Puerto Rico
4452 x 1378
5.9
4.65 min
93 min
Sierra Nevada
3750 x 2672
9.5
19.22 min
18.2 hours
Hawaii
6784 x 4369
28.2
22.35 min
Lower New England
9148 x 8509
77.8
114 min
Panama
11283 x 10862
122.5
3.5 hr
killed after 6 days
< 1% done
Outline
 A GIS application: TerraFlow
 I/O-efficient graph algorithms
• Problems and results
• Graph separation
• Topological sort on planar DAGs
 Future directions/open questions
I/O-Efficient Graph Algorithms

Input: G = (V,E)
•
Assume edge-list representation of stored on disk
Adj(v1) Adj(v2) Adj(v3) …
G

Basic problems:
•
•
•
•
BFS, DFS, CC, SSSP, MST
Standard internal memory algorithms for these problems use O(E) I/Os
Hard in external memory!
Lower bound: Ω(min{V, sort(V)})
(practically Ω(sort(V))
BFS and DFS
DFS(u)
 Mark u
 For every v in Adj(u)
• If v not marked DFS(v)
Internal memory: O(V+E)
External memory:
 one I/O per vertex to load adjacency list  O(V ) I/Os
 one I/O per edge to check if v is marked  O(E) I/Os
 O(V+E)= O(E) I/Os
I/O-Efficient Graph Algorithms

Problems:
1. Random (unstructured) accesses to the adjacency lists of
vertices as they are visited  Ω(V) I/Os
2. Need to check if v has been already visited and/or read its
key  Ω(E) I/Os
•
•
o(E) algorithm: solve (2)
o(V) algorithm: solve (1) and (2)
Upper Bounds
Sparse Graphs
 Sparse graphs E=O(V)
• CC, MST: O(sort(V)) if graph stays sparse under edge contraction
• Undirected BFS:
O(
V
)
B
• Undirected SSSP:
O(
V
B
• Undirected DFS:
O(V)
Directed BFS, DFS, SSSP
O(sort(V)) ? open
log
W
) O(sort(V)) ? open
w
o(V) ? open
 O(sort(N)) BFS, SSSP, (DFS) on special classes of sparse graphs
• Planar
• Outerplanar, grid, bounded-treewidth
Planar Undirected Graphs
 BFS, DFS, SSSP: O(sort(N)) I/Os
• O(sort(N)) I/O-efficient reductions [ABT’00, AMTZ’01]
• Separators can be computed in O(sort(N)) I/Os [MZ’02]
DFS
O(sort(N)) I/Os [AMTZ’01]
O(sort(N)) I/Os [ABT’00]
BFS
separators
O(sort(N)) I/Os [ABT’00]
SSSP
Planar graph separation:
R-division
 A partition of a planar graph using a set S of separator vertices into
. ( N ) subgraphs (clusters) Gi of at most R vertices each such that:
O
R
N
) separators vertices in
R
• There are
total
• There is no edge between a vertex in Gi and a vertex in Gj  i  j
• Each cluster is adjacent to O( R) separator vertices
O(
R
R
R
R
R
R
R
R
R
R-divisions and
Planar Graph Algorithms
 R-divisions [Frederickson’87]
dynamic graph algorithms [GI’91,KS’93], faster SP algorithms
[HKRS’97], SP data structures
 In external memory choose R = B2
• O(N/B) separator vertices
• O(N/B2) clusters of O(B2) vertices each and O(B) boundary vertices
• Can be computed in O(sort(N)) I/Os [MZ’02]
 B2-division
SSSP, BFS, DFS, topological sort, APSP, diameter,
SP data structures,..
Directed Ear Decomposition (DED)
 A directed ear decomposition of a graph G is a partition of G into
simple directed paths P0, P1, …, Pk such that:
• P0 is a simple cycle
• endpoints of each Pi i>0 are in lower-indexed paths Pj, Pl, j,l<i
• internal vertices of each Pi i>0 are not in any Pj j<i
 G has a directed ear decomposition if and only if it is strongly
connected (exists directed cycle containing each pair of vertices u,v).
 Planar DED: O(sort(N)) I/Os [ATZ’03]
Planar Topological Sort using DED
 Theorem [KK’79]: The directed dual of a planar DAG is strongly
connected and therefore has a directed ear decomposition.
 Idea:
• Place vertices to the left of P0 before vertices to the right
• Sort two sets recursively
 Used in PRAM topological sort algorithm [KK93,K93]
 PRAM simulation  O(sort(N)log N) I/Os
 Improved to O(sort(N)) by defining and utilizing ordered ear
decomposition tree [ATZ’03]
O(sort(N)) Topological Sort
using B2-division
 Construct a substitute graph GR using B2-division
• edge from v to u on boundary of Gi
if exists path from v to u in Gi
GR has O(N/B2)· O(B2)=O(N) edges, O(N/B) vertices
 Topologically sort GR (separator vertices in G):
• Store in-degree of each vertex in list L
• Maintain list of in-degree zero vertices
• Repeatedly:
• Number an in-degree zero vertex v
• Consider all edges (v,u) and
decrement in-degree of u in L
 Incorporate vertices in G in the topological order of GR
B2
v
O(sort(N)) Topological Sort
using B2-division
 Problem:
• Not clear how to incorporate
removed vertices from G in
topological order of separator
vertices (GR)
5
A
1
E
D
F
C
43
B
2
 Solution (assuming only one in-degree zero vertex s for simplicity):
• Longest-path-from-s order is a topological order
• Longest paths to removed vertices
B2
t
locally computable from longest-paths
to boundary vertices
s
O(sort(N)) Topological Sort
using B2-division
1. Compute a B2-division of G
2. Construct substitute graph GR using
• Weight of edge between v and u on boundary of Gi equal to
length of longest path from v to u in Gi
3. Compute longest path to each vertex in GR (same as in G):
• Maintain list L of longest paths seen to each vertex
• Repeatedly:
v
• Obtain longest path for next
vertex v in topological order
• Consider all edges (v,u) and
update longest path to u
4. Find longest path to vertices inside clusters
O(sort(N)) Topological Sort
using B2-division
Analysis
 Access to adjacency list of each vertex takes O(N/B) I/Os
 But..need dist(s,u)
for all u in Adj(v)

..need indegree(u) for all u in Adj(v)
 Keep list LS={dist(s,u), for any u in S}
• For each vertex v read from LS the current distances of adjacent
vertices
 O(N) edges => O(N) accesses to LS  O(N) I/Os
R-division
 Boundary vertices Bnd(Gi) of Gi
• The separator vertices adjacent to Gi
 Boundary set
• Maximal subset of separator vertices that are adjacent to the same clusters
 Lemma [Frederickson’87]:
N
• R-division of a planar graph of bounded degree has O ( ) boundary sets.
R
O(sort(N)) Topological Sort
using B2-division
Longest paths on GR
Idea: use boundary sets
Store LS so that vertices in the same boundary set are consecutive
• There are O(N/B2) boundary sets
• Vertices in same boundary set have same O(B) neighbors in GR
assuming G has bounded degree
• Each boundary set is accessed once by each neighbor in GR
• Each boundary set has size O(B)
 O(N/B2) x O(B) = O(N/B) I/Os
Open Problems/Future Directions
 I/O-Efficient Graph Algorithms
•
•
•
•
O(sort(N)) DFS on planar digraphs
Planar DAGs: can a DFS-tree be computed using topological order?
Improved algorithms for general digraphs
Simple and feasible O(sort(N)) algorithms for planar graphs and in
particular for triangulations
 I/O-Efficient GIS
• Processing LIDAR data
• Point to grid conversion, point to TIN conversion, terrain
simplification, Delaunay triangulation…
• TINs
• Practical algorithms on triangulations
• Flow modeling on TINs
End
Upper Bounds
 General undirected graphs
VB
O
(sort(
E
)
log
log
)
• CC, MST:
E
[MR’99, ABT’01]
• BFS:
VE
O(
 sort( E ))
B
• SSSP:
O(
• DFS:
E
O ((V  ) log V sort( E )) [KS’96]
B
VE
W
log
 sort( E ))
B
w
[MM’02]
[MZ’03]
 General directed graphs
• BFS, DFS, SSSP: O((V  E ) log V sort( E )) [BVWB’00]
B
Topological sort
Planar DAGs
Summary and Open Problems
 If the B2-division is given
• Topological sort can be computed in O(scan(N)) I/Os
• Extends to BFS and SSSP
 Simplified O(scan(N)) algorithms for planar DAGs
B2-division scan(N) SSSP
BFS
scan(N)
?
Topological sort
?
DFS
Planar SSSP
B2
1. Compute a B2-division of G
t
s
2. Construct a substitute graph GR on the separator vertices such that it
preserves SP in G between any u,v in S
• replace each subgraph Gi with a complete graph on Bnd(Gi)
• for any u, v on Bnd(Gi), the weight of edge (u,v) is δGi(u,v)
GR has O(N/B2)· O(B2)=O(N) edges and O(N/B) vertices
3. Compute SSSP on GR
4. Compute SSSP to vertices inside clusters
Planar SSSP
SSSP on GR with O(N/B) vertices and O(N) edges
 Dijkstra’s algorithm with I/O-efficient p-queue
• Access to adjacency list of each vertex takes O(N/B) I/Os
• O(N) Insert/Delete/DeleteMin in O(sort(N)) I/Os [A95]
• But..need dist(s,u) for all u in Adj(v)
v
 Keep list LS={dist(s,u), for any u in S}
• For each vertex v read from LS the current distances of adjacent vertices
 O(N) edges => O(N) accesses to LS  O(N) I/Os
Planar SSSP
SSSP on GR
Idea: use boundary sets
Store LS so that vertices in the same boundary set are consecutive
• There are O(N/B2) boundary sets
• Vertices in same boundary set have same O(B) neighbors in GR
assuming G has bounded degree
• Each boundary set is accessed once by each neighbor in GR
• Each boundary set has size O(B)
 O(N/B2) x O(B) = O(N/B) I/Os
Planar APSP



Straightforward bound: O(N sort(N)) = O(sort(N2))
Improved to optimal O(scan(N2))
Idea: compute SP from all vertices in a cluster while cluster is in memory
 For each cluster Gi
 For any α in Bnd(Gi) compute SSSP(α) in GR
 For each cluster Gj
 load in memory Gj, Bnd(Gj) and δ(Bnd(Gi), Bnd(Gj))
 compute the shortest paths between all vertices in Gi and Gj
d(u,v)=min{δGj (u,α) + δGR(α, β) + δGi(β,v) | α in Bnd(Gi), β in Bnd(Gj)}
α
u
Gi

v
β
Gj
 write the output
 O(N/B2) clusters  O(sort(N2)/B) [compute] + O(scan(N2)) [output]
Diameter: O(sort(N2)/B)
General AP-BFS
The APSP idea (compute SP from all vertices of a cluster while the
cluster is in main memory) can be generalized to other algorithms
which use clustering, like the BFS algorithm [MM’02] on general
undirected graphs.
Theorem:
• AP-BFS of a general undirected graph and its unweighted
diameter can be computed in O(V sort(E)) I/Os.
Note:
• general undirected BFS is O(sort(E)) amortized over V vertices
Planar DFS
Idea: Partition the faces of G into levels around a source face containing s
and grow DFS level-by-level
• Levels can be obtained from BFS in dual graph
• Structure of levels is simple (bicomps are cycles)
• Rooting/Attaching: use that a spanning tree is a DFS-tree if and
only if it has no cross edges
A DFS-tree of a planar graph can be computed in O(sort(N)) I/Os
s
s
s
H0
H1
H2
Planar Graphs
 Shortest paths
• generalize to digraphs: compute B2-division on the underlying graph
• BFS, SSSP in O(sort(N))
• APSP (transitive closure) in O(scan(N2))
• diameter in O(sort(N2)/B)
 DFS
• Undirected
• O(sort(N)) using BFS in the dual
[O(sort(N) log N) direct algorithm using cycle separators]
• Directed
• The planar undirected DFS algorithms do not extend to digraphs
• O(sort(N)) DFS? open